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ABSTRACT 

Compressive loads can cause local buckling in composite laminates that have a 
near surface delamination. This buckling causes load redistribution and secondary 
loads, which in turn cause interlaminer stresses and delamination growth. The 
goal of this research effort was to enhance the understanding of this instability- 
related delamination growth in laminates containing either an embedded or an 
edge delamination. 

There were three primary tasks: 1) development of a geometrically nonlin- 
ear finite element analysis named NONLIN3D; 2) performance of a parametric 
analytical study to determine the effects of strain, delamination shape, and de- 
lamination size on the distribution of the strain-energy release rate components 
along the delamination front; and 3) performance of a combined experimental and 
analytical study of instability-related delamination growth (IRDG). Two material 
systems (AS4/PEEK and IM7/8551-7) and two stacking sequences (0/90/90/0) 6 
and (90/0/0/90)6 were examined. The laminates were fabricated with Kapton 
inserts between the fourth and fifth plies from the top surface to give an initial 
delamination. 

* The analysis predicted a large variation of Gj and G // along the delamination 
front. The Gjjj component was always small. The location of maximum Gj 
and Gjj depended on the delamination shape and applied strain. In general, 
the strain-energy release rates were small except in a small region. Hence, 
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delamination growth was expected to occur over only a small portion of the 
delamination front. Experiments corroborated this prediction. The laminate 
stacking sequence had a large effect on the shape of the deformed region, the 
direction of delamination growth, and the strain at which delamination growth 
occurred. These effects were predicted by the analysis. The Gj component 
appeared to govern initial delamination growth in the IM7/8551-7 laminates. 
Matrix ply cracking generally accompanied delamination growth. In some cases 
fiber microbuckling also occurred shortly after delamination growth occurred. 
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Labels for nodes used in G calculation 
Semi-axes of elliptical delamination in x- and in- 
directions, respectively 
Constitutive coefficients 
Young’s moduli for orthotropic material 
Nodal forces 

Mode I, mode II, mode III, and total strain-energy 
release rates 

Shear moduli for orthotropic material 
Thickness of sublaminate 
Thickness of base laminate 
Jacobian matrix 

Coefficients in tangential stiffness matrix 

Shape function for node m 

Number of Gaussian quadrature points 

Number of shape functions 

Transverse load 

Nodal displacements 

Perimeter coordinate 

Displacements in x-, y-, and z-directions 

Strain energy 

Volume 

Transverse displacement in center of plate 
Width of finite element model. 

Also used to designate Gaussian quadrature 
weighting coefficients. 
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NOMENCLATURE, concluded 


x,y,z 
A a 

Si 

e J 

Vl2 > ^ 23 » ^13 

n 


Rectangular Cartesian coordinates 
Increment in delamination length 
Local coordinates in parent element, i= 
Strains 

Poisson’s ratios for orthotropic material 

Total potential energy 

Stresses 


Subscripts and Superscripts 

i, j =1,6 except as noted 

a, 0 = 1, number of nodal displacements 

6 =1, number of Gaussian quadrature points 



Chapter 1 
INTRODUCTION 


In recent years there has been an an increasing interest in the application of 
advanced composite materials to structures. These materials exhibit high specific 
strength and stiffness. Through the choice of fiber/matrix systems and lamination 
stacking sequences, there are great opportunities to tailor the material stiffness and 
thermal expansion coefficients in various directions to meet specific design needs. 
Some are very resistant to corrosive environments. There are also advantages 
in fabrication of certain structures, in which the “part count” can be drastically 
reduced by using composites. The list of desirable characteristics is long. 

Unfortunately, there are also potential problems associated with the use of 
composite materials. One of the primary problems is lack of experience. New 
material systems are being introduced on a regular basis. Even most of the “older” 
material systems are only a few years old. Application of composites to strength 
critical structures has been very limited. 

There are also new modes of failure, such as delamination, fiber breakage, and 
intralaminar cracking. These materials may be very strong in certain directions, 
but they can also be surprisingly weak (compared to metals) in other directions. 
The use of high strength, low strain to failure fibers severely reduces plastic 
deformation. Whereas a metallic structure can experience local yielding to reduce 
stress concentrations, a composite structure can do little readjusting to reduce 
locally high stresses without local failure. 

To expedite the process of safe application of composites to strength critical 
structures, much research has been aimed at developing appropriate stress anal- 
yses. This is not to say that metal structures do not require sophisticated stress 
analyses. But, for the same level of geometric and loading complexity, composites 
are much more demanding of the analyst than are metals. 

There are basically two types of stress analyses which must be performed. The 
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first consists of analysing a problem in which the failure is well documented and 
understood. The analyst’s task is simply to quantify the magnitude of certain 
parameters which are known to control the failure mode. For many situations it 
is not clear what governs failure. Since it is not possible to model all aspects of 
the configuration, approximations must be made. This is the second type of stress 
analysis problem. 

The problem addressed by this thesis is of the second type. Failure of a 
laminated composite material under compression loads is a very complicated 
process. Besides all of the failure modes commonly observed under tensile loads, 
there are additional failure modes related to instability. There is the possibility of 
fiber microbuckling, lamina buckling, and global buckling of the entire laminate. 
Usually, the concern is that some initial damage or defect will precipitate the 
operation of one of these instability-related mechanisms. The initial defect may 
be quite simple: perhaps a disbond between two lamina. Initial damage can be due 
to high stresses during service or impact damage. Impact damage is a particularly 
messy situation for the analyst. It is impossible to model all of the fiber breaks, 
delaminations, and lamina cracks. The task is to make approximations which make 
the analysis tractable without losing any of the essential elements of the problem. 
Learning to identify the essence of a complicated problem involves dissection of 
the original problem into independent, less complicated problems that exhibit one 
or two failure mechanisms. These less complicated problems are then carefully 
analyzed to determine the basic behavior. Hopefully, certain mechanisms can be 
identified as unlikely or at least unimportant. Then potential interactions of the 
mechanisms would be considered. 

Presently, there is no consensus as to what are the critical mechanisms of 
compression failure. Likely, different modes are critical under different situations. 
Also, more understanding of the basic mechanisms of compression failure is 
needed. This thesis will examine the failure mode instablity-related delamination 
growth; that is, delamination growth which is caused by localized buckling of 
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a delaminated group of lamina. This buckling causes load redistribution and 
secondary loads, which in turn cause interlaminar stresses and delamination 
growth. The next several sections will survey the literature and describe the 
scope of this investigation. 

1.1 Survey of the State of the Art 

Most of the research on instability-related delamination growth (IRDG) has 
used one of two configurations: the through-width delamination or the embedded 
delamination. These are shown in Fig. 1.1 . When a critical compressive strain 
level is reached, the delaminated region buckles. This causes interlaminar stresses 
along the delamination front, possibly leading to delamination growth. The term 
“sublaminate” will be used to refer to the buckled group of plies. The term “base 
laminate” will be used to refer to the unbuckled group of plies. These regions are 
labeled in Fig. 1.1. 

The literature on instability-related delamination growth differs in the choice 
of configuration, the type of stress analysis, and the method of characterizing the 
magnitude of ; he delamination front stresses. Also, some papers are primarily 
experimental and others concentrate on just the analysis. The literature survey in 
the next two sections will be organized according to the configuration considered. 

1.1.1 Through-Width Delamination 

The primary motivation for considering the through-width delamination is 
that it is less complicated than the embedded delamination. Hence, it provides a 
convenient vehicle for checking various ideas about modeling. 

Kachanov was perhaps the first to analyse the through-width delamination 
(ref. 1). He developed an approximate nonlinear beam-column analysis for the 
case of thin-film buckling. The change in strain-energy in the thin film due to 
delamination was compared with the rupture energy. If sufficient strain-energy 
was released, the delamination was assumed to grow. Kachanov did not consider 
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Through-width Delamination 




Fig. 1.1 Two configurations which exhibit instability-related delamination growth. 
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the individual components of strain-energy release rate. Whitcomb (ref. 2,3) 
developed a geometrically-nonlinear finite element analysis for calculating the 
mode I and mode II strain-energy release rates (G/ and G//). Linear finite 
element analysis was later combined with nonlinear beam theory to obtain results 
for many configurations (ref. 4, 5). This hybrid analysis drastically reduced the 
computational effort required to perform a parametric study. Chai, et al (ref. 6) 
developed a beam-column analysis for calculating total strain-energy release rate 
( Gj' = G / + G//). Gt was calculated based on differentiation of the total strain 
energy. This analysis could handle the case of global bending combined with local 
buckling. The analysis could also handle multiple delaminations, although no 
provision was made to prevent interpenetration of the different groups of buckled 
lamina. Wang (ref. 7) performed a finite element analysis of composites with 
multiple delaminations. Constraints were included to prevent interpenetration 
of adjacent lamina groups. However, only the bifurcation buckling problem was 
studied. Since there was no postbuckling analysis, strain-energy release rates could 
not be calculated. Ashizawa (ref. 8) also published a beam-column analysis. He 
calculated Gj based only on the moment at the end of the delamination. Simitses, 
et al used plate analysis to study how local buckling of a delaminated group of 
lamina affects global stability (ref. 9). In ref. 10 Sallam and Simitses presented 
a plate analysis which can be used to calculate the effects of coupling between 
bending and stretching on Gt • 

Experimental measurements of instability-related delamination growth have 
been published for both fatigue (ref. 3, 4, 11, 12) and static (ref. 11, 12, 13) loads. 
No unusual experimental techniques were used, so individual references will not 
be discussed. 

1.1.2 Embedded Delamination 

Kachanov (ref. 1) was also perhaps the first to analyse the embedded 
delamination. He presented a thin-film analysis for a circular delamination in 
a plate subjected to uniform radial loads. Prediction of delamination growth was 


7 


based on the change in strain energy in the buckled thin film versus the rupture 
energy. The individual modes of strain-energy release rate were not considered. 
Konishi and Johnston (ref. 14) assumed a trigonometric form for the transverse 
displacements of a rectangular delamination with an initial imperfection. The 
moments along the delamination boundary were estimated by differentiation of 
the assumed displacement function. Results from this analysis were used to design 
specimens for a study of delamination growth under fatigue loads. However, no 
analytical results were presented. Also, experimental results were not compared 
with the analytical model. Chai (ref. 15) developed a Rayleigh-Ritz analysis for 
the postbuckling of an embedded elliptical delamination. The buckled group of 
lamina was assumed to be very thin compared to the rest of the laminate. The 
entire laminate was assumed to be isotropic. The total strain-energy release was 
obtained by differentiation of the strain-energy with respect to the lengths of the 
axes of the ellipse. This procedure gave an average measure of the strain-energy 
release rate along the boundary. The analysis was later modified to be able to 
analyse orthotropic laminates (ref. 16). Webster (ref. 17) developed a Rayleigh- 
Ritz analysis for determining bifurcation buckling of a circular delamination. He 
used Ashton and Whitney’s technique (ref. 18) to account for bending-extension 
coupling of unsymmetric sublaminates. Shivakumar and Whitcomb also presented 
a bifurcation buckling analysis (ref. 19). They used both finite element and 
Rayleigh-Ritz plate analyses. Orthotropic, elliptical sublaminates with various 
orientations relative to the load direction were considered. Since postbuckling 
was not considered, there was not the possibility of a strain-energy release rate 
calculation in either ref. 17 or 19. Whitcomb and Shivakumar (ref. 20) presented 
a technique for calculating the distribution of strain-energy release rate along the 
delamination front from plate analysis results. Although this technique is general, 
only results for isotropic laminates were presented. Yin (ref. 21) developed a plate 
analysis for calculating total strain-energy release rate for axisymmetric loading 
of an isotropic plate with a circular delamination. Fei and Yin (ref. 22) developed 
a similar analysis for axisymmetric bending of an isotropic plate. 
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All of the preceding analyses are for quasi-static loading. Bottega (ref. 23) 
developed a dynamic plate analysis for estimating the effects of dynamic loading 
on delamination growth. Only axisymmetric configurations were considered. 

A variety of experimental programs have observed growth of embedded de- 
laminations. Konishi and Johnston (ref. 14) used Kapton implants to create 
delaminations of a known size. They performed both static and fatigue tests. 
Rhodes, et al (ref. 24) observed that impact of compression panels could lead to 
delamination, followed by local buckling and delamination growth. Byers (ref. 25) 
and Porter (ref. 26) studied the growth of delaminations which originated either 
due to impact or an implant. Ramkumar (ref. 11) also performed tests on lam- 
inates with implants. Chai et al (ref. 27) combined high-speed photography and 
shadow moire to study dynamic growth of delaminations due to combined impact 
and compression. 

1.2 Scope of Investigation 

The literature survey showed that there has been no detailed analysis of any 
truly three-dimensional ( 3D ) configuration which exhibits instability-related 
delamination growth. Even the approximate analyses that have been performed 
have examined only a very limited range of parameters. None of the analyses 
have been capable of calculating the individual modes of strain-energy release 
rate. The goal of this investigation is to enhance the understanding of instability- 
related delamination growth through detailed stress analysis supplemented by 
experiments. 

There are two facets of the analytical study: development of a suitable stress 
analysis and using it to calculate strain-energy release rates. 

Geometrically nonlinear 3D analysis is required to perform the detailed analy- 
sis. In particular, finite element analysis will be used. Such analyses are inherently 
expensive. However, geometrically nonlinear 3D analysis has become more prac- 
tical since the introduction of supercomputers, which have processing capabilities 
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in the neighborhood of 200 million floating point operations per second. To fully 
exploit the capabilities of a supercomputer requires tailoring of one’s program to 
fit the architecture of the computer. Programs written for ordinary mainframe 
computers are not generally suitable for supercomputer usage. Hence, one of the 
tasks in this study was to develop a finite element program which exploits a su- 
percomputer architecture. This program is named NONLIN3D. The development 
and verification of NONLIN3D will be discussed. 

There are several aspects to the use of NONLIN3D. Even with the use of 
supercomputers and special purpose programs, efficient modeling must be used. 
This means avoiding excessive mesh refinement and using nonlinear analysis only 
where it is necessary. Also, tasks like mesh generation can become intractable 
unless approached properly. These aspects of the use of NONLIN3D will be 
discussed. 

NONLIN3D was used to perform a parametric study of two 3D configurations 
which exhibit instability-related delamination growth. In particular, the embedded 
delamination (Fig. 1.1) and the edge delamination (Fig. 1.2) will be examined. 
A few results will be presented for the through-width delamination, which is 
basically two-dimensional. The parameters studied include strain, delamination 
size, and delamination shape. Also, the effect of stacking sequence was studied. 
Delamination growth behavior was predicted based on the calculated strain-energy 
release rates. 

Another use of NONLIN3D was to determine how accurate plate analysis is for 
calculating total strain-energy release rates. Plate analysis is potentially attractive 
because it is inherently much less expensive than 3D analysis. Results from the 
literature were used for the plate analysis. (No plate analysis was performed as 
part of this study.) 

Some preliminary experiments were performed on a configuration which does 
not buckle, but does exhibit geometric nonlinearity. This configuration was used to 
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Fig. 1.2 Laminate with postbuckled edge delamination. 


check the degree of agreement one might expect between analysis and experiment 
for relatively simple configurations. 

Experiments were performed to determine whether instability-related delam- 
ination growth could be predicted using strain-energy release rate parameters. 
Laminates with embedded and edge delaminations were tested. The deformation 
and growth of the buckled region were monitored during the tests. Also, X-rays 
and light microscopy were used to determine the extent of delamination and other 
types of damage. Predicted and observed behaviors are compared. 

The following chapters will begin with a discussion of NONLIN3D. Then the 
results of a parametric analysis of homogeneous quasi-isotroptc lamtnates with 
either an embedded or edge delamination will be presented. Then the experimental 
procedure will be discussed. Finally the results of a combined analytical and 
experimental parametric study will be presented. 
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Chapter 2 

FINITE ELEMENT ANALYSIS 


This chapter describes the theoretical aspects of the finite element analysis 
program N0NLIN3D. This program was developed as part of this thesis effort to 
perform the stress analyses presented in later chapters. There were two primary 
reasons for developing a new program rather than using a commercially available 
one. The first reason is related to the size of the computational task. Nonlinear 
3D analyses generally require large computer resources, both in terms of memory 
and number of computations. Supercomputers, such as the CDC VPS-32 and 
the CRAY-2, are well suited for the task. To exploit the power of such machines 
requires special programming techniques. Unless a program is written specifically 
for a supercomputer, the program will usually not perform well. NONLIN3D was 
written to exploit the capabilities of the VPS-32, which was the supercomputer 
available for this work. The second reason for developing a new program was to 
permit tailoring of the source code to suit the needs of this particular research 
effort. This tailoring involved all aspects of the analysis, including ease of input 
and output, simple techniques for specifying boundary conditions, and automated 
strain-energy release rate calculations. 

The following topics will be covered in this chapter: 

2.1 Governing Nonlinear Equations 

2.2 Finite Element Approximation 

2.3 Numerical Integration 

2.4 Vectorized Implementation 

2.5 Eigenvalue Analysis of Stiffness Matrix 

2.6 Substructuring 
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2.7 Contact Analysis 


2.8 Strain-Energy Release Rate Calculation 

2.9 Material Properties 

2.l Governing Nonlinear Equations 

This subsection discusses the derivation of the equilibrium equations and the 
expressions for the internally generated nodal forces and the tangent stiffness 
matrix. Also, the strain-displacement relations are discussed. 

The total potential energy IT is given by 

n = i J CijS^dV -F a q a (2.1.1) 

where the integral term is the strain energy and the second term is the potential 
energy of the applied loads. The terms C,y and are terms in the constitutive 
matrix and the strains, respectively. The terms F a and q a are the generalized 
forces and displacements, respectively. The adjective “generalized” is used to 
indicate that F a and q a need not be nodal forces and displacements in the usual 
sense. For example, in traditional Rayleigh-Ritz analyses, the q a are simply 
unknown coefficients in the series expansion for the displacements. However, 
in this discussion the F a and q a will always refer to the nodal forces and 
displacements in the z-, y-, and z-directions. The system is assumed to be 
conservative; hence, the equilibrium state is obtained by minimizing IT, which is 
accomplished by setting the first partial derivatives with respect to the unknowns 
equal to zero. 


dU 

dq a 


I 


Cij^i 


dej 

dq f* 


dV - F a -0 


( 2 . 1 . 2 ) 


Equation 2.1.2 is nonlinear because of the nonlinear strain-displacement rela- 
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tions. The integral in eqn. 2.1.2 gives the magnitude of the internally generated 
nodal forces corresponding to the current displacements. Until a converged so- 
lution is obtained, the internally generated forces do not equal the externally 
applied forces. The differences in the forces, referred to as residuals, equal 
The Newton-Raphson procedure that was used to solve eqn. 2.1.2 requires the 
partial derivatives of the residuals with respect to the unknowns. These partial 
derivatives of the residuals are the coefficients in the tangential stiffness matrix K 
and are given by eqn. 2.1.3. 


K a(3 = 


a 2 n 

dq a dqP 



dcj 

dq@ 


dV + 


/ 




d 2 e,- 

dq a dqP 


dV 


(2.1.3) 


The first integral gives the terms for the sum of the linear and large dis- 
placement matrices. The second integral gives the terms in the geometric stiffness 
matrix. If the strains are equal to zero, the matrix obtained using the first integral 
and the nonlinear strain-displacement relations is identical to the matrix obtained 
by simply updating the coordinates and using the linear strain-displacement rela- 
tions. If the strains are not equal to zero, there is a difference. This was verified 
numerically for 2D elements. 

A Lagrangian formulation was used in NONLIN3D. For infinitesimal strain the 
nonlinear strain displacement relations are (ref. 28) 


ei = Ux + 1/2 (u x u x + V X V X + W X W X ) 

£2 = Vy + l/2(u y Uy + VyVy + WyWy) 

C3 = Wz + 1/2 {u z u z + v 2 v z + w z w z ) ^ ^ 

£4 = Uy d* V X (tixUy "H V X Vy W X Wy) 

£5 = V 2 + Wy + (u Z U y + V Z Vy + W Z Wy) 

£6 = u z + w x + ( U Z U I + v t v x + WzWx) 
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where u, v, and w are displacements in the x — , y — , and z-directions, and the 
subscripts x, y, and z indicate partial differentiation (e.g. u y = ^). The nodal 
values of u, v, and w are the unknowns referred to earlier as q a . Note that 6:4, 
£■5, and £6 are engineering shear strains. Since a Lagrangian formulation is used, 
the strains are based on the original configuration. For example, £1 is the axial 
strain of a line which was originally (i.e., before deformation) parallel to the x- 
axis. Although this line could be oriented parallel to the y-axis after deformation, 
the axial strain is still £1 (not e 2 )■ 

The material coefficients, C t} , were assumed to couple normal and shear strains 
in the xy plane only; hence, C^j = 0 for j = 1,2, 3, 4, and 6 and C$j = 0 for j = 
1 through 5. This corresponds to an orthotropic material which has one material 
axis parallel to the z direction. In expanded form, the stresses are given by 


'C11 

c 12 

C13 

C14 

0 

0 - 


■ffl' 

c 12 

C‘ 22 

c 23 

c 24 

0 

0 


£ 2 


c 23 

<?33 

c 34 

0 

0 


£3 

c 14 

C 24 

C 34 

C44 

0 

0 


£4 

0 

0 

0 

0 

<?5S 

0 


£5 

. 0 

0 

0 

0 

0 

c 66 - 


.£-6- 


2.2 Finite Element Approximation 

Equations 2.1.1 through 2.1.5 are general equations which can be used with any 
procedure based on minimization of total potential energy. In the finite element 
method, the body is divided into subregions referred to as elements. Within 
an element, the displacements u, v, and w are approximated by interpolation 
functions N m and nodal values of the displacements, u m , v m , and w m . 


u = N m u m 


v = N m v m where m = 1, NS 

[NS = number of shape functions) 

w = N m w m 


( 2 . 2 . 1 ) 
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The shape functions for 8- and 20-node elements are given in the appendix. 
Fig. 2.2.1 shows schematics of the two elements. 

Since an isoparametric formulation was used, the element geometry is 
approximated using the same interpolation functions and the nodal coordi- 
nates, z m , y m , and z m . 


x — N m x m 


y = N m y m 


( 2 . 2 . 2 ) 


2 = N m z m 


Calculation of the tangential stiffness matrix (eqn. 2.1.3) requires the first 
and second partial derivatives of the strains (eqn. 2.1.4) with respect to the 
nodal displacements u m , v m , and w m . Since the nodal values are independent of 
the coordinates, the partial derivatives with respect to coordinate directions take 
the form u m . Also, the nodal displacements are independent. Hence, 

= £mn and, Qftpr = 0, where 6 mn is the Kronecker delta. Equation 2.2.3 
gives the expressions for the derivatives of the strains. 
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( 2 . 2 . 4 ) 


The matrix B is introduced to define the terms B tr , which will be used later 
to conveniently refer to the derivatives in eqn. 2 . 2 . 3 . Most of the second partial 
derivatives are zero. The non-zero ones are given in eqn. 2 . 2 . 5 . Note that the 
terms are the same for derivatives with respect to u, v, and w. 
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d 2 ci 

du^du* 


d 2 €2 

du m du* 


d 2 ej 

du m du** 


d 2 64 

du^du* 


d 2 £p 

5u m 5u^ 


d 2 e Q 

du^du^ 



- _ 

dv m &v n 

dw m dw n 


d 2 €2 _ d 2 e 2 _ 

dv m dv n ~ dw m dw n ~ 


d 2 e s 

_ d 2 e s _ 

dv m dv n 

dw m dtv n 

d 2 C4 

~ 3 2 e 4 - 

<?v m dv n 

dw m dxv n 

<9 2 €T5 

— d2e j - — 

dv m dv n 

dw m dw n 

d 2 ea 

— _ 

5v m dv n 

dw m Bw n 


N?N? 


K N y 


N?N? 


N™Ny + N™N? 


N™N? + N™N£ 


N™N? + N™N2 


= a7 in 


— a mn 
~ °2 


= a 


a 


mn 

3 

mn 


= a; 


mn 


= a? n 


(2.2.5) 


2.3 Numerical Integration 


The integrals in eqns. 2.1.2 and 2.1.3 were evaluated using numerical integra- 
tion. The integrations become summations as shown in eqns. 2.3.1 and 2.3.2: 


NG 

(2-3.1) 
(2.3.2) 

where 

NG = number of quadrature points 

| J| = determinant of the Jacobian (required because integrations 
are performed using local coordinate system) 




dq c 


6-1 


dq 


NG 


K ‘* = 


NG 


d^£i 


W = Gaussian quadrature weighting coefficients 

and the superscript 6 indicates that the term was evaluated at quadrature 
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point 6. 


The first summation in eqn. 2.3.2 can be expressed in a symmetric form. 
This form is like that in ref. 29, but the motivation is different. Herein, the 
symmetry is the valuable characteristic which will be exploited in the vector 
implementation of the finite element analysis. The material coefficients are 
replaced by C xj = Q is Q ja , where Q ta = 0 for i < s. The terms Q ta are obtained 
using Cholesky decomposition. Cholesky decomposition is always possible since 
the constitutive matrix C x} is symmetric and positive definite. The weighting 
coefficients W and the determinant of the Jacobian | J\ are positive, so the square 
root of \J\W is real (only the positive square root need be considered). Integration 
schemes which involve negative weighting coefficients are not considered herein. 
Hence, the first summation in eqn. 2.3.2 can be written as 


NG 

K? 

8=1 


(2.3.3) 


and T aa is defined to be 

T’ a = Qi,§^{\J\Wy 5 s = l,NSTR 

(2.3.4) 

where NSTR — number of strains 

The product is now that of the transpose of a matrix and itself. 

For the 8-node element, the terms related to normal strains were evaluated at 
8 points (i.e. a 2 x 2 x 2 integration scheme). To improve the performance of 
the 8-node element in modeling bending deformation, terms related to er 5 and er 6 
were evaluated at the centroid only (ref. 3, 30, 31). The other shear strain terms 
related to e 4 were evaluated with full integration. Both 2x2x2 and 3x3x3 
integration schemes were considered for the 20-node element. There was not much 
difference in the behavior of the element for the two integration schemes, so the 
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2x2x2 integration scheme was used for the results presented herein. 

2.4 Vectorized Implementation 

The three-dimensional finite element program NONLIN3D developed for use 
in this study was designed to exploit the vector processing capabilities of a 
supercomputer, the CDC VPS32. The CDC VPS32 is closely related to the CDC 
CYBER 205. In the context of computation on the VPS32, a vector is simply 
a list of numbers; it is not a vector in the usual mathematical sense, wherein 
a vector has magnitude and direction. The architecture of the VPS32 is such 
that long vectors can be manipulated (i.e., multiplied, added, etc.) more than 
an order of magnitude faster them the individual terms could be manipulated 
separately. Ref. (32,33) gives some quantitative information on the processing 
speed. Processing speed is often expressed in terms of millions of floating point 
operations per second (MFLOPS). Actually, the speed is not the same for different 
operations. For example, multiplication is performed considerably faster than 
division. The following rates are applicable for multiplication and addition. For 
scalar operations, the processing speed is approximately 3-5 million MFLOPS. For 
very long vectors, the speed is 100-200 MFLOPS. Short vectors are not processed 
nearly so fast. However, a vector length of 110 to 160 will result in a rate of 50-100 
MFLOPS. Obviously, it is highly advantageous to manipulate long vectors rather 
than scalars or short vectors. 

To vectorize a program simply means to implement the mathematical algo- 
rithms in such a way that long lists of numbers are manipulated rather than 
individual numbers. There is usually no unique vectorization of a particular task. 
There is considerable room for ingenuity. In fact, successful vectorization of a task 
generally requires discarding of procedures familiar to the scalar programmer (for 
example, those procedures in familiar textbooks) and taking a fresh look at what 
has to be accomplished. 

There are three subroutines which perform essentially all of the computation 
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intensive work in NONLIN3D. These are BOBL3D, ESTIF3D, and the equation 
solver. The subroutine BOBL3D calculates strains, stresses, and element forces. 
The subroutine ESTIF3D calculates element stiffness matrices. The equation 
solver will not be discussed, since vectorized versions of banded (ref. 34) and profile 
(a vectorized version of that in ref. 35) equation solvers are already available. 

Reference 36 describes a vectorized procedure for calculating element stiffness 
matrices for linear high-order elements, which require a large number of quadrature 
points. The primary technique in ref. 36 was to manipulate the values of a 
parameter at all the quadrature points simultaneously, rather than individually. 
This resulted in vectors of length equal to NG (the number of quadrature points) 
and, at times, NG * NS, where NS = number of shape functions. At one stage of 
the stiffness matrix calculation in ref. 36 even larger vectors were manipulated, 
but the longer vectors were obtained at the expense of significant replication of 
vectors and fairly complicated logic. Reference 37 describes a technique which 
results in long vector lengths for low order elements, but it has the same drawbacks 
mentioned above for the technique in ref. 36. Also, the technique in ref. 37 has 
very large memory requirements. Techniques were developed in the current study 
which result in fairly long vector lengths without much replication, complexity, 
or memory requirements. Also many of the techniques used can be extended to 
process multiple elements at the same time, thereby increasing vector lengths. 
Simultaneous multiple element processing was not attempted for two reasons: 
l) time constraints and 2) the expected reduction in total computation costs did 
not appear to justify further optimization of the routines. 

To expedite the discussion of BOBL3D and ESTIF3D, special notation will be 
introduced. Also, simple vector programming methods are discussed. Vectors are 
given abbreviated names. A bar under a name indicates that it is the name of a 
vector. When vectors are separated by the symbol “ * ”, vector multiplication is 
implied. For example, 
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(2.4.1) 


A * B = C 


is equivalent to A x B t = C x , with no sum on i. The vectors may also be stacked, 
for example, 


'a' 

B 

* 

'D' 

E 


A * D 

B*E 

C 


F 


C*F 


(2.4.2) 


In all cases, vector multiplication refers to term by term multiplication without 
any summation. A dot “ ■ ” will be used to denote an inner product. A VPS32 
special function is used to perform the inner products as a single vector operation. 

It is very convenient to manipulate vectors just like scalars. For example, 
suppose we need the determinant of ten 2x2 matrices. In vector form, we would 
write the matrices as follows 


A 

C 

In eqn. 2.4.3 all 10 values of a coefficient are grouped in a single vector. Hence, 
the length of each vector in the matrix is 10. The 10 determinants are obtained 
using two vector multiplications and one vector subtraction 


B 

D 


(2.4.3) 


the ten determinants = A*D_ — C_*B_ (2-4.4) 

Using scalar methods there would have been 20 scalar multiplications and 10 
scalar subtractions. 

Often it is desirable to perform replication of a scalar or vector in order to 
reduce the number of subsequent vector operations. For example, suppose we 
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need 20 copies of the vector a stacked contiguously in the vector D_. Consider D_ 
to consist of 20 subvectors d\ each of which have the same length as a. 

The first step is to assign values to d 1 ; that is, set d 1 = a. Then assign d 2 = a. 


then assign 

[d 3 


Id 1 } 

d 4 


d 2 


rd 5 i 


■d 1 ! 


d 6 


d 2 

Then assign 

d 7 

= 

d 3 


.d 8 J 


_d 4 J 


Then assign 


r d 9 1 


rd 1 1 

d 10 


d 2 

d 11 


d 3 

d 12 


rf 4 

d 13 

— 

d 5 

d 14 


d 6 

d 15 


d 7 

Ld 16 , 


Ld 8 J 


Finally assign 


rd 17 i 


rdH 

00 


d 2 

d 19 

II 

d 3 

Ld 20 J 


Ld 4 J 


By using this technique the number of vector assignments is reduced from 20 
to just 6. 

The next two subsections describe the two routines B0BL3D and ESTIF3D. 

2.4.1 Subroutine B0BL3D 

Figure 2.4. 1.1 shows a flowchart for B0BL3D. There are 8 primary tasks 
performed by this routine. The first task is to calculate local derivatives of 
the global coordinates at each of the quadrature points. The shape functions 
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Fig. 2.4. 1.1 Flowchart for subroutine BOBL3D. 




are expressed in terms of local coordinates S r (see appendix), but derivatives 
with respect to global coordinates are required. The first step towards obtaining 
these derivatives is to calculate the local derivatives of the global coordinates, 
Slf-’ an> ^ Since an isoparametric formulation is used, the geometry is 

approximated with the same shape functions as the displacements. For example, 
x = N m x m where x m = the x coordinate of the m tri node. Hence, 


dx d{N m x m ) _ dN m 
36^ ~ Wr ~ dS T 


(2.4.1.1) 


There are NS shape functions to be differentiated at NG quadrature points. 
For a single element type, the are invariant, so they are calculated only once 

and stored in DL in the following order. 



dN 


N 1 I 


3? 7 


N 2 

DL = 

dN 

where N_ = 

N 3 


dN 

db$ 


.N» s . 


and N_ m = the m t ^ shape function evaluated at each of the NG quadrature points. 
Note that N m has a length of NG. (The length of DL is 3 * NS * NG.) 

The nodal coordinates for an element are stored in vectors X, Y, and T±. Each 
nodal coordinate is replicated NG times. The ordering is as follows for the vector 
X. 
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x = 



(2.4. 1.3) 


Note that the vector X has a length of NS * NG. The vectors Y and Z axe 
ordered the same way. The vectors X, Y, and Z are stacked and replicated once to 
obtain a long vector XXYYZZ . This stacking and replication reduces the number 
of vector multiplications required in the next step. The organization of XXYYZZ 
is given by eqn. 2.4. 1.4. 


XXYYZZ = 



(2.4. 1.4) 


The vector XXYYZZ is of length 6 * NS * NG. 

The derivatives of the global coordinates are obtained in two steps. First, three 
vector multiplications between DL and parts of XXYYZZ are performed and the 
results are stacked in as shown in eqn. 2. 4. 1.5. 
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(2.4.1.5) 


The length of D is 9 * NS * NG. 

As indicated in eqn. 2. 4. 1.1, there are NS terms to be summed for each of 
the nine partial derivatives. After the summations, there will be NG values of 
each partial derivative. The vector D in eqn. 2.4. 1.5 is organized such that there 
are nine subvectors, each of length NS * NG. Recursive addition is used on each 
subvector to sum the appropriate terms. The vector product of (which is part 
of DL) and X will be used to illustrate the recursive addition. The number of 
nodes, NS, will be assumed to be 8. This product constitutes the first subvector 
in eqn. 2.4. 1.5. The vector product of and X creates a vector which consists 
of the subvectors d t (see eqn. 2. 4. 1.6). 



X 

DL* 

Y 


Z 

'y' 

DL* 

z 


X 


z z z 

DL * 

X 


Y 



[d l l 

d 2 

d 3 

d 4 

d 5 

/ 

d 7 

Ld 8 J 


(2.4. 1.6) 


where d m — the contribution from shape function m. The length of each d m is 
NG. 

For example, d 2 = * x 2 evaluated at each of the NG quadrature points. 

Now three recursive additions are performed to add up the d m . 
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(2.4. 1.7) 


-<*r 


r<*r 


r^i 

d .2 

da 


—2 

da 

+ 

d& 

in 

-^4- 


-di- 


■is . 


then 


'i r 



+ 

is 

.—2 


.—2 

kJ 


(2.4. 1.8) 


finally, = <L\ + d^ 

The vector d± now contains evaluated at NG quadrature points. The 
other eight derivatives ( etc.) are obtained in a similar fashion. For 

example, the derivative is obtained by replacing with 62 and X with Y in 
eqn. 2. 4. 1.6. 


The second task in B0BL3D is the calculation of the determinant and inverse of 
the Jacobian at each of the quadrature points. The Jacobian consists of derivatives 
of the global coordinates with respect to the local coordinates. The scalar form is 
shown in eqn. 2. 4. 1.9. 



dx 

dy 

dz 

d 6 \ 


d 6 \ 

dx 

062 


dz 

dS 2 

dx 

dy 

dz 

56 3 

do$ 

Ml 


(2.4. 1.9) 


In the vector implementation, the entries in J are vectors (the vectors which 
were calculated in task 1). For example, J_n is a vector consisting of evaluated 
at NG points. The inversion and determinant are performed explicitly in terms 
of these vectors. The result is an inverse which has entries which are vectors of 
length NG and a determinant vector of length NG. 
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Next the global derivatives of the shape functions are calculated at each 
quadrature point (task 3, Fig. 2.4. 1.1). These derivatives are calculated by 
multiplying the local derivatives by the inverse of the Jacobian. The scalar form 
is shown in eqn. 2.4.1.10. 


dN m , 
dx 


-u ii 

IJ 12 

U\3 ' 


r dNfn n 

~MT 

dN m 

3y 



IJ 21 

IJ 22 

IJ 23 


dN m 

dN™ 

dz 


IJ 31 

IJ 32 

IJ 33 


dN m 


where the matrix IJ is the inverse of J. 


(2.4.1.10) 


In the vector implementation, the entries in eqn. 2.4.1.10 are replaced by 
vectors of length NG. For example, the NG values of are obtained as follows 


dN_ 2 

dx 


Lin 


dN 2 
~ds T 


+ Hl2 * 


dN 2 

dh 


+ Hl3 * 


dN 2 
dd 3 


(2.4.1.11) 


where / = the i, j vector term in the inverse Jacobian matrix, calculated in 
task 2, and N} is defined in eqn. 2.4. 1.2. 

The nine global derivatives of the displacements ( ••••etc.) with respect 

to the global coordinates can now be calculated (task 4). A typical term is: 


du _ _ dN m 

~d~x~ Ux ~ dx 


* u 


(2.4.1.12) 


For convenience, a subscript is used to indicate the variable of differentiation. For 
example, u x = and u z = 

This is exactly the same form as in eqn. 2.4.1. 1. The local derivative 
has been replaced by the global derivative and the nodal coordinates, x m , 
have been replaced by nodal displacements, u m . Consequently, the same vector 
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procedure that was used earlier to calculate local derivatives of coordinates is now 
used to calculate global derivatives of displacements. A vector DG is assembled 
from the global derivatives of the shape functions. 


DG 


f N r l 


N y 


L Nj 


(2.4.1.13) 


where N is defined in eqn. 2. 4. 1.2. The subscript indicates the variable of 
differentiation. Note that the length of the vector DG is 3 * NG * NS. 

Equation 2.4.1.13 is simply a “global derivative” version of eqn. 2. 4. 1.2. Also, 
a displacement vector uuwww is assembled. 


(2.4.1.14) 


The organization of u, y, and w are identical to x, y, and z respectively. For 
example, in eqn. 2. 4. 1.3 replace x*, x 2 , ... by u 1 , u 2 ,... Now the vector procedure 
used earlier can be used by simply replacing DL by DG and xxwzz by uuwww . 
Two copies of the shape function derivatives, N x , Ny, and N z , are stored in 
DSXYZ (eqn. 2.4.1.15). 


uuwww - 


u 

V 

w_ 

u 

V 

. w . 


DSXYZ = 


DG 

DG 


(2.4.1.15) 


The vector DSXYZ is of length of 6 * NG * NS. 

Task 5 is to calculate strains. The strains are calculated using eqn. 2.1.4, 
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except that the scalar variables are replaced by vectors of length NG. Thus, the 
value of a strain component is calculated at all integration points simultaneously. 
For example, 


£x = M* + 2 * ~ x + - 1 * ~ x + ^ ^ 


(2.4.1.16) 


This is possible because the global derivatives of the displacements (calculated 
in task 4) are organized such that for each displacement derivative, all NG values 
are contiguous in memory. The strains are stored in the vector e in the following 
order. 


£ = 


£i 

£2 

£3 

£4 

£5 

-£ 6 - 


(2.4.1.17) 


where the length of each e, is NG. 

Task 6 is to calculate the partial derivatives of the strains with respect to 
the displacements. The scalar form of these derivatives is shown in eqn. 2.2.3. 
These involve derivatives of the shape functions (for example, N and products 
of the derivatives of the shape functions and displacements (for example, u z N m ). 
There are 27 products. To reduce the number of vector operations in forming the 
products, stacks of vectors are manipulated, as described next. 

The displacement derivatives calculated in task 4 are replicated (for a total of 
NS copies of each derivative) and stacked in the vector DD . 
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DD = 


’EH* ' 

gVy 

Rtu g 

Rv z 

RW y 

Enz 

R™X 

Ru 

Lflv'J 


Ru x = 


u. 


u. 


NS copies 


similarly for other subvectors in DD. 


(2.4.1.18) 


Note that the displacement derivative vectors u x , v y , etc. are replicated NS 
times. Hence, the vector DD is of length 9 * NS * NG. This is in preparation for 
multiplication with the shape function derivatives, since, for example, the length of 
u x is NG but the length of N y is NG * NS. To expedite the discussion of the use of 
these replicated vectors, the use of the prefix “R” will be used to denote replication. 
The required products are calculated by performing vector multiplications of parts 
of DD with parts of DSXYZ . These products are stored in a vector BLV . 
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BLV 


* ElLx " 


[Kx] 

RVy 


Ky 

Rw r 

£ 

Nj 

Rv r 


Nj. 

Rw y 


Ny 

- Elkz - 


Ljv'J 

Env i 


r^i 

Rw z 


Ky 

Enx 

£ 

Kz 

RWy 


N x 

Ru 7 


Ky 

-EHU- 


LjvJ 

' Rwj ‘ 


r^xi 

R\U 


Ky 

RWy 

* 

Kz 

RU; 


Kz 

Ryu 





L AT, J 


RV y * Nj 


Rw r 


'Kz 

EUy 

★ 

Ky 

_Rv z . 


UJ 


Rux * N_ y 


Ry y 

* 

'Kz' 

. Eiu . 


[Ky J 


Ru„ * N 7 


Rri z * N_ x 


(2.4.1.19) 


Note that only 9 vector products are required to form the 27*NS + NG products 

in BLV. 
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Now it is a simple matter to form a vector of the 18 strain derivatives by adding 
vectors stored in DG and BLV. 


■B n- 


Kx + «*JV* 

M 21 


UyKy 

Msi 


HzKz 

M 41 


Ky + OxKy + « y Kx 

Msi 


HzKy + HyKx 

Msi 


N z + H t N_x + u r^r 

Ml2 


VxKx 

—22 


Ky+HyKy 

M 32 


v*Kz 

Ma2 

— 

Kx + V-yKt + V t N y 

Ms 2 


Kz + Vy Kz + V-zKy 

M 62 


VxKz + VxKx 

—13 


u^Nx 

B.23 


—y—y 

K 33 


K z + w z N z 

Maz 


WyN_X + W T .Vy 

Ms 3 


Ky + W-zKy + WyKz 

Ms 3 


Kx + aixKi + uuK x 


where the B.^ are vector versions of the B{j shown in eqn. 2.2.4. The length of 
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each is NS*NG. Since the B,y are formed using vector operations, the organization 
of each B^y in terms of the NS shape function contributions and NG quadrature 
point values is the same as for N, described in eqn. 2.4. 1.2. For example, 


fin* 


B\i 


B h 


B? 


11 


mX s 


and the length of each subvector B™ is NG (2.4.1.21) 


Task 7 is to calculate the product of the stresses, determinant of the Jacobian, 
[Jj, and weighting coefficient W. The first step is to assemble a vector containing 
copies of | J| *W evaluated at the NG quadrature points. 


DJW 

DJW 


DDJJWW = 


DJW 


NSTR replicates of DJW (2.4.1.22) 


where DJW = | J\ * W evaluated at NG quadrature points and NSTR = number 
of strains (NSTR = 6 for 3D analysis) 

The vectors of the strains £, are stacked in memory, (see eqn. 2.4.1.17) so 
that only a single vector multiplication of the stacked strains and DDJJWW is 
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required to obtain the scaled strains e . 


■ir 


'£1' 

£2 


£2 

£3 

= DDJJWW * 

£3 

£4 


£4 

£5 


£5 

L£e J 


L£6-l 


(2.4.1.23) 


The scaled stresses are obtained by linear combination of the scaled strains. 


d* = Cijij (2.4.1.24) 

This is analogous to the unsealed, scalar stresses in eqn. 2.1.5. For efficiency 
later, the d, are stacked contiguously in d in a form which is exactly analogous to 
the stacking of e t in eqn. 2.4.1.17. 

The final task in BOBL3D is to calculate the element forces R a corresponding 
to the current stress state. The element forces are calculated using a vectorized 
version of the summation in eqn. 2.3.1. For clarity in the following discussion, the 
displacements q a are separated into u, v, and w, the displacements in the x-, y-, 
and z-directions. The vector form of the strain derivatives in eqn. 2.3.1 are in the 
large vector B (see eqn. 2.4.1.20). 

The vector implementation of the summation in eqn. 2.3.1 will be illustrated 
by describing the calculation of the x-direction forces next. The first step in 
calculating the x-direction forces is to evaluate the products in eqn. 2.3.1 at all 
NG quadrature points and for all NS nodes. Recall that values of the triple product 
(< 7 *|J|W)* are already available in the vector d . At each quadrature point each 
stress must be multiplied by NS values related to each derivative of the strains. To 
reduce the number of vector multiplications, the weighted stresses are replicated 
and stored as follows. 
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where Rd. is a replicated form of dj 


Rd = 


Mi' 

M2 



For example , 



NS copies of dj 



(2.4.1.25) 


The length of Rd is 6 * NS * NG 


Now all the products for the x-direction can 
multiplication (eqn. 2.4.1.26). 


be formed with a single vector 


r/ii 


rfini 


-*r 

L 


8.21 


£2 

Lz 


Hn 

* 

£3 

L 


B 41 


£4 

U 


B 51 


£s 

UJ 


-Sj&\ - 


-£e- 


(2.4.1.26) 


Now recursive addition is used to sum the /^, as follows. 


’Li 


h 


[£.1 

—2 


L 2 

IL\ 

+ 

Lie J 


(2.4.1.27) 


then /j = /j + f_2 + Lz' 
The content of /. is now 


L = 


NG terms to be summed to obtain F J 
NG terms to be summed to obtain F~ 


. NG terms to be summed to obtain F* s J 


(2.4.1.28) 
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where F™ = x-direction force at node m. 

The summations required to obtain the F™ from the terms in / are performed 
using a VPS32 special function which is optimized for summing contiguous 
numbers in storage locations. The y- and z-direction nodal forces are obtained 
in a similar manner. 

2.4.2 Subroutine ESTIF3D 

The routine ESTIF3D is used to calculate the tangential stiffness matrix, which 
consists of the large displacement matrix, K i and the geometric (or initial stress) 
matrix (eqn. 2.3.2). The large displacement matrix involves derivatives of the 
strains (which are in B^ } ), the Cholesky factors Q< a for the constitutive matrix 
(Cij = QisQjg), and the square root of the products From eqn. 2.3.3, we 

see that for each term {Ki ) a ^ , there are NG inner products to be evaluated. For 
example, 


NG 

K}* = Y^{T sl T aZ ) 6 s = l,NSTR 

6=1 


(2.4.2.1) 


The vector T is formed which contains the T am at the NG quadrature points 
and is organized such that all NG inner products for a particular stiffness matrix 
term can be performed at once. The first step in forming the vector T is to perform 
a vector version of the product Q{ a * The strain derivatives are available in 
the vectors J9 tj . Since the Cholesky factor terms Q{ a are scalars, the resultant 
product vectors P_ t} are simply linear combinations of the vectors B tJ . 


Pi, = Qt&j (2. 4.2. 2) 

The organization of the product vectors P,y is the same as for B tJ (in terms 
of contributions for NS nodes and NG quadrature points). Also, the P,, are 
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stacked contiguously in memory as P just like the in 5 (see eqn. 2.4.1.20). 
Now the vector is reordered so that terms which will summed are contiguous. The 
reordering is performed using a special VPS32 function referred to as Q8VGATHR. 
Equation 2. 4. 2. 3 gives the order of T , which is the reordered version of P. 



similarly for 
shape functions 
2 through NS 
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In T, the P™ vectors are placed according to the shape function with which the 
term is associated. The same is the case for other P™ vectors. Recall that in P T J , 
all NS P™ vectors were stacked together. The subvectors T™, T™, and T™ are 
introduced here for convenience later. These subvectors are stored contiguously in 
eqn. 2. 4. 2. 3 for convenience later in determining pointers for the inner products. 
The length of each subvector is 6 * NG. The subscripts x, y, and z indicate that 
these subvectors contain terms related to the x-, y-, and z-directions, respectively. 

The terms in T must be multiplied by the square root of the determinants 
and weights (i.e. y/\J\W terms) to obtain T. To reduce the number of 
vector multiplications, the square root of DJW (which has a length of NG) was 
calculated and then replicated to obtain the vector SDJW which has a length of 
3 * NS * NSTR * NG. Then only a single vector multiplication is required to form 
the vector T . 


(2.4. 2.4) 

The terms in Ki are simply the inner products involving the subvectors PJ 1 , 
TJ 1 , and V?. 

For example, the term in K \ related to the first shape function and the x- 
direction and the fourth shape function and the y-direction is 


f = T* SDJW 


lltl ( 2 . 4 . 2 . 5 ) 

The inner products (indicated by the a • ”) were performed by a VPS-32 special 
function. 

The geometric stiffness matrix involves much less manipulation than Ki. The 
terms C X} oAJ\W are already available in the vector &. Equation 2.2.5 gives 
the scalar form of the second partial derivatives of the strains. Note that the 
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derivatives are independent of the direction of the displacements variable. Hence, 
the three nonzero stiffness terms associated with each coordinate direction for each 
combination of shape functions are the same. These derivatives consist of products 
of the derivatives of the shape functions with respect to the global coordinates, 
which are stored in the vector DSXYZ . For a particular shape function, the 
NG values of each second partial derivative is calculated by performing vector 
multiplications and additions with the appropriate subvectors (of length NG) in 
DSXYZ . The second partial derivatives are stored contiguously in memory, so 
each stiffness term is calculated with a single inner product of the stress vector a 
and the second partial derivatives. 


k mn = 


a™ 


—2 


n mn 


,mn 


_mn 


_mn 

“6 


■& 


( 2 . 4 . 2 . 6 ) 


, mn _ d 2 g« _ 

where flj du m du n dv m dv n dw m dw n ‘ 

The k mn are assembled in the element stiffness matrix as follows. 


rfc 11 

0 

0 

k 12 

0 

0 ..." 

0 

k n 

0 

0 

Jfc 12 

0 ... 

0 

0 

k n 

0 

0 

k 12 ... 

k 12 

0 

0 

k 22 

0 

0 ... 

0 

Jb 12 

0 

0 

k 22 

0 ... 

0 

0 

1 k 12 

0 

0 

Jk 22 ... 

; 

! 


• 

• 

» ' « 


( 2 . 4 . 2 . 7 ) 
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2.4.3 Execution Times 


This section discusses the results of several CPU “time trials” which quantify 
the advantages of using vector algorithms on the VPS32. The following discussion 
is for the 20-node element and a 2 x 2 x 2 integration scheme, except as noted. 

One of the programs available in the Fatigue and Fracture Branch at NASA 
Langley Research Center used a non-vectorized algorithm to calculate the linear 
element stiffness matrix for the 20-node element. This algorithm took .33 second 
of CPU to calculate the linear element stiffness matrix for one element. A vector 
version of this program took only .014 second of CPU for the same calculation. 
The program NONLIN3D was tailored for geometrically nonlinear analysis. No 
provision is provided (as of this writing) for calculating just the linear stiffness 
matrix. So a comparable timing is not available for NONLIN3D. However, the 
total time for NONLIN3D to obtain l) the tangential stiffness matrix (which 
includes both the large displacement stiffness matrix and the geometric stiffness 
matrix), 2) the element strains and stresses, and 3) the element forces was only 
.015 second per element. Since these three calculations involve many more terms 
than the linear stiffness matrix, the algorithms developed for NONLIN3D can be 
considered to be very well vectorized. 

Timings were also obtained for NONLIN3D for a 3 x 3 x 3 scheme. For the 
three tasks listed in the previous paragraph, the total time was .0215 CPU second 
per element. This is only 40 percent more time than for the 2x2x2 scheme. For 
a scalar code, the time per element would have increased by more than a factor of 
three, since there are 27 quadrature points for the 3x3x3 scheme instead of the 8 
for the 2x2 scheme. The reason that the time increased by only 40 percent is that 
the more refined integration scheme has longer vector lengths, but the number of 
vector operations is unchanged. 
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2.5 Eigenvalue Analysis of the Element Stiffness Matrices 


The finite element stiffness matrix should permit rigid body translation and 
rotation of the element without inducing any strain- energy. Further, there should 
be strain-energy whenever the element is deformed. For displacement formulated 
finite elements which are integrated with exact numerical integration, there is 
usually no problem. However, at times it is desirable to use a lower order 
of numerical integration in order to improve the performance of an element in 
modeling bending deformation (ref. 3, 30, 31) or simply to reduce the number of 
computations required to form the element stiffness matrix. Whenever a reduced 
order of numerical integration is used, there is the possibility of introducing zero- 
energy deformation modes; that is, the element can deform in certain modes 
without the expenditure of any work. 

A convenient method of detecting zero-energy deformation modes is to perform 
an eigenvalue analysis on the element stiffness matrix. The eigenvalue analysis is 
described next. The following discussion is based in large part on ref. 38. 

The eigenvalue problem has the general form 


K6 = A 6 


(2.5.1) 


The solution of eqn. 2.5.1 yields “n” eigenvalues and eigenvectors, where n = 
the number of degrees of freedom (DOF) in the element. For example, the 20-node 
3D element has 60 DOF, so n = 60. For a 3D element there should be 6 rigid 
body modes ... 3 translations and 3 rotations. An excess of 6 rigid body modes 
indicates the presence of spurious zero-energy deformation modes. Less than 6 
rigid body modes indicates that strain is associated even with rigid body motion. 
The eigenvalue corresponding to each rigid body mode should equal zero. This is 
apparent when eqn. 2.5.1 is pre-multiplied by the transpose of 6. 
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S t K6 = 6 t X6 


(2.5.2) 


The left side of the equation is simply twice the strain energy in the element 
when it undergoes the displacements 6. For rigid body motion the energy should 
be zero. The right hand side of the equation is zero only if all the displacements 
are zero (which is the trivial solution ) or if A = 0. 

Early check-out runs of NONLIN3D used the reduced integration 8-node 
element to analyse a very simple configuration. The stiffness matrix was found 
to be singular, which indicted that the zero-energy deformation modes might be 
causing problems. Eigenvalue analyses of the element stiffness matrix showed 
that there is one zero-energy deformation mode for each shear strain treated 
with reduced integration. The cause for this can be seen readily by considering 
a set of displacements which causes a shear strain S5 which varies linearly in 
the x-direction, does not vary with y or z, and is zero at the centroid. If 
only the centroidal value of the strain is sampled, the strain will be determined 
to be zero and there will be no energy associated with the deformation. For 
many problems the boundary conditions are such that the spurious zero-energy 
deformation modes do not cause any major problems. Section 4.2 discusses the 
results of an attempt to use the reduced integration 8-node element for analysis 
of a postbuckled sublaminate. 

2.6 Substructuring 

The program NONLIN3D was designed to perform analysis by substructures. 
A brief description of the substructuring technique is given here. More details 
can be found in ref. 39. In addition to reducing computer memory requirements, 
substructuring allows the structure to be modeled as a combination of linear and 
nonlinear components. For the configurations studied herein (figure 1.1 and 1.2), 
linear analysis is appropriate everywhere except the majority of the postbuckled 
region. By substructuring into linear and nonlinear regions, expensive iterative 
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solution is needed for only a fraction of the equations. 

In this study, two substructures were used: one linear and one nonlinear. The 
substructuring procedure begins by obtaining a reduced stiffness matrix and load 
vector for the linear region. The reduced stiffness matrix can be treated as the 
stiffness matrix for just another type of element. Because of the large number of 
nodes, this element will be referred to as a superelement. The stiffness matrix 
and load vector are “reduced” in the sense that only the nodes shared by the two 
substructures (the interface nodes) are included. 

The technique used for calculating the reduced stiffness coefficients utilized the 
formal definition of a stiffness coefficient: a stiffness coefficient is related to the 
restraint forces required to maintain unit displacement at one degree of freedom 
(DOF) and zero displacement at the remainder of the element DOF. Suppose there 
are to be n DOF in the superelement. These n DOF are restrained. One of these 
DOF is specified to have a unit displacement (and there are no other loads) and 
the governing equations are solved. The restraint forces at all n DOF constitute 
one column of the reduced stiffness matrix. This procedure is repeated for all n 
DOF. The reduced load vector is obtained in a similar fashion. All n DOF in 
the superelement are still restrained. However, now the specified loads for the 
linear region are applied. The reduced load vector is equal to the negative of the 
restraint forces at the restrained nodes. 

Once the superelement stiffness matrix and load vector are calculated, the 
analysis proceeds to the nonlinear substructure. Whenever the nonlinear stiffness 
matrix and load vector are formed, the interaction with the linear substructure 
is included by simply adding the superelement stiffness matrix and load vector. 
When the internally generated nodal forces are calculated to determine residuals, 
the contribution of the linear substructure consists of the product of the superele- 
ment stiffness matrix and the superelement nodal displacements. 

For the configuration analyzed, the delamination front is within the linear 
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substructure. Hence, further work is required even after obtaining a converged 
solution for the nonlinear substructure. After obtaining a converged solution, 
the displacements for the interface nodes are known. These displacements fully 
account for the effect of the nonlinear substructure on the linear substructure. 
That is, the displacements in the linear substructure can be determined as 
though there was no other substructure, except that the magnitudes of the 
displacements at the interface nodes are specified. To reduce the computer 
resource requirements, it is usually advantageous to obtain multiple solutions 
for the nonlinear substructure and then obtain multiple solutions for the linear 
substructure. 

2.7 Contact Analysis 

For certain combinations of delamination size and strain level, closure occurs 
over part of the delamination front. Unless constraints are imposed, the delam- 
ination faces will overlap in the analysis. When contact constraints are added, 
there are two types of geometric nonlinearity: that due to significant rotations 
and that due to the unknown contact area. Since the contact area affects the 
global response and vice versa, the iterative loop includes two inner iterative loops 
for obtaining either the effects of significant rotation or the contact area. Because 
of the selected substructuring, there are two substructures in which contact can 
occur. But, as mentioned in section 2.6, one of the substructures was treated with 
linear analysis. Strictly speaking, there is not a linear substructure now. However, 
an “engineering” (as opposed to a purist) approach was taken in implementing the 
contact analysis. If only a small fraction of the original interpenetration is per- 
mitted, then most of the effect of imposing full contact constraints will be seen. 
Also, because of the nature of the problem, large interpenetration cannot occur 
in the “linear” substructure if large interpenetration is prevented in the nonlinear 
substructure. Hence, the approach taken was to perform nonlinear analysis on the 
buckled region only. 

There was also an approximation made in how the constraints were imposed. 
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Fig. 2.7.1 will be used to explain this approximation. The figure shows a disk 
resting on a base. Constraints are imposed to prevent the disk from overlapping the 
base. Ideally the constraints would be imposed such that the constraint force on 
the disk is provided by the node it contacts on the base. A free-body diagram would 
look like Fig. 2.7.1b. In NONLIN3D the constraint force on the disk is supplied 
externally, so the free-body diagram looks like that in Fig. 2.7.1c. Obviously there 
is some error due to using an external force, but the amount of error and its 
importance depends on the individual problem. 

Unfortunately, for the postbuckled sublaminate, the computational effort 
would have been too large to perform a more rigorous analysis to determine the 
amount of error. However, since the base laminate has a w=0 boundary condition 
on the plane z=0, the effect of applying an external force rather than reacting the 
contact force on the base laminate is probably small. 

The flow chart in Fig. 2.7.2 outlines the procedure used to perform the 
contact analysis. In the flowchart the term nonlinear solution refers to solving 
the governing equilibrium equations (eqns. 2.1.2) with the current set of contact 
constraints. The first step is to check for w displacements which would cause 
overlap. Since the w displacements of the base laminate are so small, it is sufficient 
to search for nodes with negative w displacements. If there are none, the current 
solution is correct. For each node with a negative w, a constraint and a load 
are imposed such that the node will have w=0. Next, a nonlinear solution is 
obtained for the current contact constraints and applied load. Then the signs of 
the contact forces are checked. Any node with a tensile constraint force must have 
the constraint released. If there are no tensile restraint forces, the current solution 
is correct. After releasing all tensile constraints, another nonlinear solution is 
obtained. Now the loop begins again with checking for negative w. 
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Any nodes with w < 0 ? 




no 


yes 


For all nodes with w < 0: 

* add constraint at node in z-direction 

* add corrective load at node to make 
w — 0. 


Iterate to obtain new nonlinear solution 
for current set of constraints 


1 

1 

Any tensile contact restraint forces ? 


yes 


Release tensile restraints 


Iterate to obtain new nonlinear solution 
for current set of constraints 


Fig. 2.7.2 Flowchart for contact analysis. 




2.8 Strain-Energy Release Rate Calculation 


The well-known virtual crack closure technique (ref. 40) served as the basis 
of the strain-energy release rate calculation. This procedure determines G/, G/y, 
and G/// from the energy required to close the delamination over a short distance, 
Aa. The closure energy involves products of delamination front nodal forces and 
relative displacements behind the delamination front. The delamination front 
nodal forces can be determined by actually closing the delamination over Aa. 
Another technique, which requires only a single solution, assumes that the current 
delamination front nodal forces are the same as they would be if the delamination 
length was reduced by A a. The single solution method was used herein. 

The strain-energy release rate calculation will be illustrated for the 20-node 
element, since this element was used for all of the parametric analysis. Figure 
2.8.1 shows a schematic of the delamination front region. The nodes of interest 
for the strain-energy release rate calculations are indicated by the filled circles. 
Because it is not appropriate to close the delamination over part of an element, 
there are four sets of nodes (indicated by the letters a, b, c, and d) which are 
used to calculate the closure energies. The relative displacements are obtained by 
subtracting the displacements at nodes a - and b ^ from the displacements at nodes 
a* and 6 t , respectively. The forces are equal to the nodal forces transmitted across 
the delamination plane at nodes c,- and d,. These forces are obtained by evaluating 
the integral / G<yc,J|i d.V for all elements which are connected to nodes c, or d* 
and whose centroids lie above the delamination plane. There are two sets of energy 
products. One of the sets of energy products consists of the relative displacements 
for nodes a, and a' t multiplied by the forces for nodes c,. The other set of energy 
products consists of the relative displacements for nodes 6 t and b ^ multiplied by 
the forces for nodes d*. The energies equal 1/2 of these products. 

Strain-energy release rate is a measure of energy per unit area. Hence, the 
energy products must be normalized by the appropriate areas. Unfortunately, 
there is not a simple exact way to determine the appropriate areas. The primary 
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complication is that the midside nodes and corner nodes are “weighted” differently 
by the assumed element shape functions. The result is that, even if the strain- 
energy release rates are actually constant along the delamination front, there would 
be much larger energy products for the midside nodes than for the corner nodes. 
For example, in figure 2.8.1 the energy products associated with nodes C2 and 
C4 would be much larger than for that associated with nodes ci and C3. An 
approximate solution to this dilemma is as follows. The strain- energy release rate 
is not calculated for locations like C2 and C4 along the delamination front. Instead, 
the energy products associated with those locations are split evenly between the 
adjacent nodes. For example, the energy associated with location C3 along the 
delamination front becomes 


~Hj — E / ■+■ , / , -f* — ~E 1 E 

a s<*s c s 2 L ®J a i c a 


a 4 a 4 c 4 J 


( 2 . 8 . 1 ) 


where E denotes the energy products associated with Gj, Gff, and Gjjj and the 
subscripts indicate the nodes involved. The area is approximated by the product 
of A a times the distance between the midside nodes on either side of the corner 
node being considered. For example, the area for node C3 is A a times the distance 
from node C2 to node c 4 . 

If the delamination front is not parallel to one of the coordinate axes, it is 
preferable to add a coordinate transformation to the procedure outlined above. In 
particular, a local coordinate system is defined for each node along the physical 
delamination front (i.e., the nodes c,). Figure 2.8.2 shows a schematic of a 
delamination plane and the global (zy) and local {x y) coordinate systems. This 
local coordinate system has one axis tangent to the delamination front, one axis 
normal to the delamination front, and one axis normal to the delamination plane. 
For all the cases considered z and z were parallel. The transformed nodal forces 
F z > , Fyi , and F g i , are 
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Fig. 2.8.2 Global and local coordinate systems and perimeter coordinate S 



( 2 . 8 . 2 ) 


F s > = F z cos 9 + F y sin 9 
F t = -Ft sin 0 + Fy cos 9 


The relative displacements are transformed similarly. The transformed forces 
and relative displacements are then used to calculate the energy products. Figure 
2.8.2 also defines the perimeter coordinate “S”, which is the distance along the 
delamination front measured from the y-axis. 

The procedure just outlined was implemented in two slightly different ways 
for the results presented here. The difference was : n the way the nodal forces 
were calculated. Initially in the study, the nonlinear strain-displacement relations 
were used to calculate the nodal forces from the nodal displacements in the linear 
region. This is inconsistent, but if the region assumed to be linear is “exactly 
linear” , it would make no difference. The results for the mesh convergence study 
were obtained using this procedure. All of the other results were obtained by 
using the linear strain-displacement relations to calculate the nodal forces from the 
nodal displacements in the linear region. Because the linear region is not exactly 
linear, there is a difference in the results obtained using the two methods. The 
configuration used for the mesh convergence study was also used in the parametric 
study, so results appear for that configuration using both methods. The second 
method is recommended by the author. 


2.9 Material Properties 


Several kinds of materials were used in this study. Some of the specimens 
involved A1 or steel. Also, two graphite composite systems were examined: 
AS4/PEEK and IM7/8551-7. The material properties used are given in Table 
1. For the two metals, the usual assumed properties were used. The assumed 
properties for the PEEK were actually some that are typical for graphite/epoxy 
(ref. 41). Since only a qualitative analysis of PEEK was performed, these 
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Table 1 Material Properties 


E n 

13.4510 

E 22 

1.02510 

E 33 

1.02510 

1/12 

.3 

^23 

.49 

"13 

.3 

Gi2 

.552510 

G 23 

.343510 

Gl3 

.552510 


* ref. 41 


IM7/8551-7** 

16.2510 Pa 
.814510 Pa 
.814510 Pa 
.22 
.22 
.22 

.648510 Pa 
.648510 Pa 
.648510 Pa 

University of Wyoming under NASA 


Generic Graphite/ Epoxy* 

Pa 
Pa 
Pa 


Pa 
Pa 
Pa 

** unpublished data generated by the 
grant NAGl-674, which began in July, 1986. 
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properties were probably sufficiently accurate. The properties for the IM7/8551-7 
were based on unpublished data generated by the University of Wyoming under 
NASA grant NAG1-674, which began in July, 1986. This grant determined the 
inplane properties, E n , E 2 2 , Gl2, and i/ l2 . Because of a lack of other data, the 
remaining 3D properties were assumed. 

In addition to the materials mentioned above, a fictitious material was used. 
For the initial parametric study, the goal was to examine the effect of strain 
level and geometric parameters on Gj, G//, and Gj/j. Material properties were 
desired which would have minimal effect on the distribution of the strain-energy 
release rates. For quasi-isotropic laminates the in-plane stiffness is independent of 
direction. But even for quasi-isotropic laminates the flexural stiffness varies with 
direction. Hence, even if the postbuckled region consisted of a quasi-isotropic 
group of plies, one would expect variations in the strain-energy release rate along 
the delamination front which are due solely to the variation in flexural stiffness. 
Also, the properties of the interface plies (i.e., those plies on either side of the 
delamination) would be expected to at least affect the percentages of G/, and 

G//, and G///. 

The simplified material properties chosen for this study are those for a 
“homogeneous quasi-isotropic laminate” throughout the entire specimen (buckled 
and unbuckled regions). These properties Cij are obtained as follows: 

(291) 

t= 1 

where ( C XJ ) l are the constitutive properties for the L tk ply in the 8-ply quasi- 
isotropic laminate (±45/0/90),. With these properties throughout, there are 
obviously no stacking sequence effects and no variation of material properties with 

orientation. 
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Chapter 3 

FINITE ELEMENT MODELING 


This chapter describes the technique used for mesh generation and discusses 
typical finite element models. 

3.1 Mesh Generation 

Figure 3.1.1 outlines the procedure used for generating most of the meshes. 
This procedure is based on the procedure in ref. 42. A two-dimensional model 
is swept through a 90° arc to generate a cylindrical 3D mesh. The outer part of 
the cylindrical mesh is then transformed to obtain a square boundary. Then an 
elliptical transformation is applied to obtain an elliptical delamination front. If 
the ellipse is longer in the y- direction than in the z-direction (i.e., b > a), the 
conformal transformation is 


x = x 


y 


= y 



b 2 — a? 
x 2 + y 2 


z = z 


(3.1.1) 


If a > b, the transformation is the same except that x and y are interchanged. 
To avoid a singularity in eqn. (3.1.1), nodes at zero radius were shifted to lie on 
an arc of very small radius, i.e. about 10“ ® m. 

The transformation in eqn. 3.1.1 maintains the orthogonality of lines which 
were orthogonal in the modified cylindrical mesh (Fig. 3.1.1c). This orthogonality 
at the delamination front simplifies the pairing of nodal forces and relative 
displacements in the strain-energy release rate calculation. 

A peculiarity of the transformation in eqn. 3.1.1 is the unusually close 
spacing of the elements close to the delamination front on the long axis of 
the ellipse (Fig. 3.1.2b). Also note what appears to be a triangular element 
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Fig. 3.1.1 Procedure for generating three-dimensional finite element models. 


c 

o 



61 


in Fig. 3.1.2c . This element has two sides which are essentially colinear. A 
modified transformation which results in more even mesh refinement is obtained 
by introducing a scale factor for the “stretching” 


y 


= y( 



b 2 — a 2 
x 2 + y 2 


F - F+ 1) 


sfx 2 4- y 2 , /— 5 5- 

where F = /or r > yx 2 + y 2 

r 


(3.1.2) 


and F = 1 for r < V* 2 + y 2 

By choosing the parameter r a little less than the radius of the delamination 
front in the cylindrical mesh, orthogonality is maintained in the neighborhood of 
the delamination front during the elliptical transformation. 

After the elliptical transformation the midside nodes are no longer at the 
middle of an element edge. Therefore, the coordinates of the midside nodes are 
recalculated as the average of the coordinates of the adjacent corner nodes. 

In section 5.3 results are presented for a square and a rectangular delaminated 
region. For the square delamination circular arcs in a mesh like 3,1.1c were 
stretched to form the side of a square. A typical mesh is shown in Fig. 3.1.3. 
For a rectangular delamination the coordinates were magnified in one direction 
to elongate the delaminated region. The transformations for the square and 
rectangular delaminations are not conformal. Hence, the orthogonality of the lines 
which are orthogonal in Fig. 3.1.1c was not preserved in the model in Fig. 3.1.3. 
This probably reduced the accuracy of the modeling. 

3.2 Finite Element Models 

In the course of this investigation the following configurations were examined: 
l) a l ami nate with an embedded delamination, 2) a laminate with an edge 
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Fig. 3.1.3 Mesh used for analysis of the square delamination. 
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delamination, and 3) a laminate with a through-width delamination. The loading 
was m-plane compression except for the transversely loaded circular delamination. 
The finite element models for these configurations are discussed in this section. 

Figure 3.2.1 shows a typical finite element model for a laminate with an 
embedded delamination. The elements are 20-node isoparametric hexahedra. 
Because of symmetry it is sufficient to model only one fourth of the specimen 
and impose the constraints u = 0 on i = 0 and v = 0 on y = 0. There is also 
a constraint w = 0 on z = 0. This constraint was imposed to remove global 
bending from the analysis. In reality, there might be global bending (particularly 
if the buckled region is thick), but the amount of global bending would depend 
on the region modeled and the boundary conditions at the external boundaries. 
Imposition of w = 0 on z = 0 simply removes overall specimen size and external 
boundary conditions as parameters to be considered in this study. The constraint 
on w represents a laminate which is well constrained globally. Of course, one could 
also view the imposition uj = 0onz = 0asan indication of symmetry about the 
z = 0 plane. This implies the presence of two delaminations. 

Along the boundary x = W, all u displacements are specified to equal We o, 
where eo is the specified compressive axial strain. To initiate transverse deflections, 
a transverse load was applied at the center of the delaminated region. After 
a converged solution w as obtained, the load was removed, and solutions were 
obtained with only compression loading. 

Figure 3.2.2 shows a typical model after division into substructures. Most of 
the postbuckled region is included in the nonlinear substructure. The distance 
between the delamination front and the beginning of the nonlinear substructure 
was L. In all cases, i was approximately equal to the sublaminate thickness h. 

The edge delamination models were similar to the embedded delamination 
models. Because the plane y = 0 is now a free surface, the constraint v = 0 on 
y = 0 used for the embedded delamination was not used for the edge delamination. 
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Fig. 3.2.1 Typical finite element model. 


Nonlinear substructure Linear substructure 
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Fig. 3.2.2 Typical division of model into substructures. 



The model for a laminate with a through-width delamination (Fig. 3.2.3) 
had exactly the same boundary conditions as the model for a laminate with an 
embedded delamination, except an additional constraint v = 0 on y = W was 
imposed to cause a plane strain response. 

A transversely loaded laminate with a circular delamination was analyzed to 
help verify NONLIN3D. The model for this configuration was like that in Fig. 3.2.1, 
except the delaminated region was circular and a transverse load (i.e a load in the 
z-direction) was applied at (x,y,z)=(0,0,H+h). 
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Chapter 4 

EVALUATION of NONLIN3D 


A variety of checks were made to verify the reliability of the finite element 
program N0NLIN3D and the models used with the program. Some of these 
checks are presented in the following subsections. 

4.1 Analysis of Transversely Loaded Plate 

A closed form solution for a circular isotropic plate subjected to a central 
point load is given in ref. 43. This solution is exact for linear deflections and 
approximate for large deflections. The equation for the central deflection w 0 is 

;.S 

(4ll) 

Finite element analyses were performed using a mesh similar to that in 
Fig.3.2.1, but with a circular debond. The thickness and radius of the debonded 
region were .4 mm and 15 mm, respectively. The Young’s modulus was 207 GPa 
and the Poisson’s ratio was .3 . 

Figure 4.1.1a compares the deflection at the center obtained from the closed 
form solution and NONLIN3D. The agreement between the two analyses is 
excellent in the linear and the initial nonlinear region. There is a 10 percent 
difference in the deflections at the highest load level considered. This is not 
surprising, since the closed form analysis is not exact for large deflections. 

The closed form solution in ref. 43 can be used to calculate the total strain 
energy release rate. Because the configuration is axisymmetric, C?r is constant 
around the boundary. The strain-energy release rate is -(dU/da)/(2-ira). The 
strain energy U can be calculated as the work done by the applied load (since all 
of the work is stored as strain-energy) . 
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Load, N 
(a) Deflection 



Load, N 

(b) Total strain-energy release rate 


Fig. 4.1.1 Analysis of a transversely loaded plate with a c 
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Equation 4.1.2 yields 


U = 2.30 


w%Eh? 

a 2 


+ .510 


w*Eh 

a 2 


(4.1.3) 


An expression for G 7 * is obtained by differentiating eqn. 4.1.3 with respect to 
a and dividing by the circumference. 


Gf = 


Ehw 2 

2xa* 


h 2 
.217 


+ 1.02 Wo 


(4.1.4) 


Figure 4.1.1b shows Gj vs. load. Results are shown for the linear closed-form 
solution (eqn. 4.1.4 without the fourth order w 0 term), the nonlinear closed-form 
solution (eqn. 4.1.4), and NONLIN3D. The results are plotted with log-log axes 
because of the wide range of the parameters. In the linear range all three solutions 
agree very well. Even after nonlinear effects become important, the nonlinear 
closed-form solution and NONLIN3D still agree very well. 

These comparisons indicate that NONLIN3D does account for geometric 
nonlinearity and that the strain-energy release rate calculation technique is valid. 

4.2 Failure of 8 -Node Element 

Initially, the 8 -node element with reduced integration was to be used for the 
stress analyses. As pointed out earlier, it was known that the element could 
exhibit spurious zero-energy modes. However, initial tests with the element seemed 
to show good performance. The critical test involved analysis of a postbuckled 
sublamanate, since this is the focus problem for this thesis. For simplicity only 
the sublaminate was modeled and constraints were applied to the lower surface 
to simulate a very stiff base laminate. As mentioned in the analysis chapter, 



postbuckling was initiated by applying both in-plane compression and transverse 
loads. After obtaining a converged solution, the transverse load was removed. 
Then the desired solutions could be obtained. For this particular test of the 8- 
node element, the transverse load was removed in steps. (In all the other cases 
presented in this thesis, the transverse load was removed all at once, not in steps.) 

Fig. 4.2.1 shows the deformed meshes for the four load cases. Converged 
solutions were obtained for all four load cases. For the first three cases, the 
deformed shape is quite smooth. All spurious oscillations are very small. When 
the last of the transverse load was removed, the deformed shape shows clearly 
that the zero energy modes are no longer subdued. A close-up of the central part 
of the plate for the case P=0 is shown in Fig. 4.2.2. These results show that 
the contribution of spurious zero-energy deformation modes depends on both the 
displacement constraints applied on the boundaries and the load system. Because 
of this behavior, it was decided that this element is not reliable. Hence, the 20- 
node element was used exclusively for the strain-energy release rate analyses. 

4.3 Check of Mesh Refinement 

Several configurations were analyzed as part of this study. Consequently, it 
was neither practical nor warranted to perform a convergence study for all cases. 
Instead, a systematic convergence study was performed for a single configuration. 
Mesh refinements f° r other configurations were selected based on the results of the 
convergence study. 

The configuration selected for the convergence study had a circular delamina- 
tion with a radius of 15 mm. The sublaminate thickness h was .4 mm and the 
base laminate thickness H was 4 mm. The overall laminate width W was 50 mm. 
The material properties were those for the homogeneous quasi-isotropic laminate 
described earlier. Figure 4.3.1 shows the extremes of refinement used for the 2D 
meshes. The elements in the coarse mesh were subdivided to obtain the refined 
mesh. Note that for the coarse mesh, only two elements are used to model most of 
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Fig. 4.2.1 Postbuckled meshes using 8-node element, (strain == -.0045) 
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Fig. 4.2.2 Closeup of postbuckled mesh (P = 0). 
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Fig. 4.3.1 Range of two-dimensional mesh refinement for convergence study. 




the buckled region. As described earlier, the 2D meshes were swept through a 90° 
arc to generate the 3D meshes. Figure 4.3.2 shows four of the meshes generated. 
As shown in the figure, the number of slices of elements was also varied. Models 
with 4 and 8 slices are shown. A 12-slice model was also used. Information on all 
of the models used in the convergence study are given in Table 2. As shown by 
the figures and the table, a fairly wide range of refinement was examined. 

Since strain-energy release rates were of primary importance, variations in G/, 
G//, and Gjjj were used to determine the adequacy of the mesh refinement. 
Figure 4.3.3 shows the distribution of G/ and G// along the delamination front 
for three strain levels and four models (models 1,3,4, and 6). Only symbols are 
shown for the crudest mesh, model 3. These meshes bracket the entire range of 
refinement in Table 2. The mode III component G//j was negligible for all cases. 
Of interest here are the differences in the results obtained using the various meshes. 
Except for model 3, which only had 4 slices of elements, the results from all of the 
models are essentially equal. Even a coarse 4-slice model gives the correct trends. 
Apparently, a fairly crude model is sufficient to calculate Gj and Gjj. Models 
with 8 slices were selected for the parametric analyses in Chapter 5, which will 
discuss the significance of the magnitude and the distribution of G/ and G//. 

The preceding convergence study was for a homogeneous quasi- isotropic lam- 
inate, which is a fictitious material. Chapter 7 presents strain-energy rates for 
actual laminates. Because a minimum of one element was required through the 
thickness of each lamina, the 2D mesh used to generate the 3D mesh had more 
nodes than the homogeneous laminate required. Unfortunately, the required com- 
putational effort prevented a systematic convergence study. In fact, to keep the 
computational effort manageable, the number of slices in the models was lim- 
ited to just four (instead of the 8 slices used for the homogeneous quasi-isotropic 
laminates). Fig. 4.3.4 shows one of the models. Based on the convergence study 
presented for a homogeneous laminate, one might estimate that the four slice mod- 
els for actual laminates underestimate the maximum Gj and Gjj by about 15-20 
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Fig. 4.3.2 Several of the three-dimensional meshes used in the convergence study. 


TABLE 2 Statistics on models used in convergence study 



QO 00 ^ 00 ^ 00 


OJ 

CO 

o 

o 

o 

OJ 

OJ 

cn 

CO 



00 



r-* 

a5 

CO 

CO 

CO 

CM 

cO 

o 

uo 

in 

in 

CM 

OJ 


co 

cd 

<Ji 





CD 

CD 


CD 


CM 

m 

in 

CO 

in 

CO 

O 

CM 

CM 

OJ 

OJ 

CO 


In 

in 


uo 


CO 

h- 

h- 

CD 

r- 

in 

r^ 



a> 


T— * 

i — 

▼— 




OJ 

in 


co 


co 

CD 


CD 

CO 

o 

CO 

CO 

o 

CM 

CM 

CO 

T— 

CO 

in 

y — 







cn 

1 — 

a> 

in 

y— 

in 



OJ 

00 


CM 


CO 


o 

o 

co 

T — 

y— 

v— 

OJ 

CO 

CD 


CD 

00 

OJ 

OJ 

CM 

OJ 

in 

CO 

00 




i— 






CO 

CO 

CO 

CO 

CO 

CO 

'O’ 

in 

1^ 

r- 


CO 

■r— 


T— 

1 — 

T~ 

in 


04 CO in CD 


>> 

a> 

> 


Q. 

CO 


Q> 


(0 

a> 

"E 

o 

c 

■D 

C 

CO 

co 

o 

c 


a> 

CO 

Oi 

■Q 

c 

CO 



78 


(a) Mode I 



(b) Mode II 

Fig. 4.3.3 Strain-energy release rates calculated using different meshes. The 
symbols indicate the results for model 3. 


Fig. 4.3.4 Typical finite element mesh used for analysis of actual laminates. 


percent. 


4.4 Determination of Allowable Residual 

NONLIN3D solves the governing nonlinear equations using a modified Newton- 
Raphson solution procedure. During each iteration the internally generated forces 
are compared with the externally applied loads. The differences are the residuals. 
If all of the residuals are identically zero, the governing equations are exactly 
satisfied. Of course, in practice exact agreement is seldom obtained. Iteration 
could continue until the algorithm’s best approximation of zero is obtained. The 
size of this “numerical” zero will depend on the computer and variable type 
specifications (i.e. single or double precision) in the program. However, negligible 
residuals are in general orders of magnitude larger than a numerical zero. 

To determine what is a negligible residual, three residual tolerances were 
considered: 1000., 1., and .0001 Newtons. A laminate with a 30 x 60 mm 
delamination was analyzed for five strain levels. The range of strains was such 
that the maximum lateral deflection for the buckled region varied from about .6 
to 2.2 times the thickness of the buckled region. The tolerance of 1000 Newtons 
gave erroneous results. The other two tolerances gave virtually identical results 
except for the lowest strain level (e x = —.001), for which there were differences of 
about 6 percent. Figure 4.4.1 shows Gf and Gjj for e x = —.001 (i.e., the worst 
case). Results for the other strains are not shown, since the differences are very 
small. 

Based on these results, a tolerance of .0001 Newton was selected for all 
the analyses. Probably a somewhat larger residual could have been tolerated. 
However, the residuals tend to decrease quite rapidly during the iterations, so 
there is little to be gained (in terms of reduced cost) by trying to specify the 
largest acceptable residual. Also, a larger tolerance might give poor results for 
some cases in which the loads are small. 
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4.5 Verification of Substructuring 


The configuration used for the convergence study was also used to verify 
the validity of the substructuring. Model 3 was analyzed with and without 
substructuring. The calculated Gj and G jj distributions from these analyses are 
shown in Fig. 4.5.1. The small difference in the results verifies the substructuring 
technique. 
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Chapter 5 

PARAMETRIC ANALYSIS OF HOMOGENEOUS QUASI-ISOTROPIC LAMINATES 

This chapter will discuss the results of a parametric study of homogeneous 
quasi-isotropic laminates containing a postbuckled embedded or edge delamina- 
tion. This fictitious homogeneous laminate was selected so that the effects of strain 
level and various geometric parameters on deformation and strain-energy release 
rate could be examined without the additional complications due to stacking se- 
quence effects. 

The following sections begin with a parametric study on the effect of strain 
level, delamination shape, and delamination size on deformation and the distribu- 
tion of G/, G//, G/// and G T along the delamination front. 

5.1 Deformation and Strain-Energy Relea se Rates For an 
Embedded Delamination 

The parameters considered were strain level, delamination shape, and delami- 
nation size. 

Fig. 5.1.1 shows lateral displacement in the middle of the delaminated region 
vs. axial strain for two circular and two elliptical delaminations. The dimensions 
of the delaminations are shown on the figure. Before buckling the lateral deflection 
is essentially zero. When the buckling load is exceeded, the displacement increases 
rapidly at first with increased strain. Then the rate of increase in displacement 
decreases. Obviously, the response is quite nonlinear. 

Fig. 5.1.2 shows plots of deformed finite element meshes for a circular and 
an elliptical delamination. The displacements have been multiplied by 10 to 
improve visualization. The deformed shape is relatively simple except near the 
delamination front. For both cases the delamination front is open near the 
intersection of the delamination front with the x=0 plane. However, for the 
circular delamination, the delamination faces actually overlap near the y=0 plane. 
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= 0 plane 
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Fig. 5.1.2 Deformed finite element meshes. Displacements scaled by 10 for (a) 
and 20 for (b). 



Even at strains smaller and larger than that shown, the circular delamination 
exhibited closing. For the elliptical delamination “small” strains result in opening 
along the entire length. At larger strains (not shown), the elliptical delamination 
also exhibited closure. Strictly speaking, constraints should be added to prevent 
overlapping of the delamination faces. However, including constraints to prevent 
overlapping further complicates an already complicated stress analysis problem. 
Consequently, no contact constraints were added for any of the results presented 
in this section. Section 5.4 presents a few results which illustrate the effect of 
including contact constraints. In the results that follow, dashed lines will be used 
for the G i and Gjj distribution curves in regions where overlap occurred. 

Figure 5.1.3 shows the Gj and G/j distributions for a circular delamination for 
five strain levels. This is the same configuration used for the convergence study 
in section 4.3. Gjjj was essentially zero for this and all other cases considered 
in this study. In general, one would not necessarily expect Gjjj to be zero. The 
strain-energy release rates are plotted in Fig. 5.1.3 using the perimeter coordinate 
5. This coordinate is zero where the delamination front meets the y-axis and is 
maximum where the delamination front meets the x-axis. Both Gj and Gjj show 
large variations along the front and are maximum at S = 0. There is overlapping 
of the delamination surfaces over a large portion of the front, as indicated in figure 
5.1.3a. Although Gj is larger than G/j for the five strain levels, the difference is 
not large; this is definitely a mixed mode situation. Since both Gj and Gjj are 
largest at S = 0, one would expect delamination growth to occur preferentially 
perpendicular to the load direction, i.e. in the y-direction. 

Since a circular delamination is expected to become elongated perpendicular 
to the load direction, a 30x60 mm elliptical delamination was analyzed. Fig. 5.1.4 
shows the distribution of G/, G//, and G t for this elliptical delamination. There is 
a large variation of both G / and G// along the front. Note that the location of the 
Gj peak shifts slightly with strain level. In contrast to the circular delamination, 
the peak values of G/ and G// occur at different locations. Also, the peak value 
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(a) Mode I 



(b) Mode II 


Fig. 5.1.3 Strain-energy release rate distributions for a circular delamination, 
(diameter = 30mm) 
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of G[j is larger than the peak value of G[ except for the case e x = —.005. The 
total strain-energy release rate (Fig. 5.1.4c) varies significantly along the front for 
the larger strains, but for the smaller strains the magnitude is almost constant. 
Depending on one’s choice of growth criterion, very different predictions of even 
the direction of growth are possible. A criterion based only on Gj would predict 
growth perpendicular or nearly perpendicular to the load direction for all of the 
strain levels. A criterion based only on Gjj would predict growth parallel to the 
load direction. For the smaller strains, a criterion based on total strain-energy 
release rate would predict almost uniform growth along the delamination front. 

Fig. 5.1.5 shows Gj and (7/j for a 60x60 mm delamination. Comparison 
of Figs. 5.1.3 and 5.1.5 show the effect of delamination size on Gj and Gjj 
for two circular delaminations. Both delaminations were subjected to the same 
strain levels. Figs. 5.1.3a and 5.1.5a show that the larger delamination is 
closed (actually over-lapping) for more of the delamination front. Also, note 
that the distribution in the overlapping region is more complicated for the larger 
delamination. This is because the strains for the larger delamination are larger 
multiples of the bifurcation buckling strain. This conclusion was verified by 
subjecting the smaller delamination to higher strains. (These results are not 
presented in this dissertation.) The larger delamination has a much larger Gj 
for the region near 5 = 0. Figs. 5.1.3b and 5.1.5b show that Gjj is also larger for 
the larger delamination near 5 = 0. Hence, one might expect unstable extension of 
the delamination once it begins to grow. However, based on the calculated strain- 
energy release rates, a circular delamination is not expected to grow self-similarly 
into a larger circular delamination. It should become elliptical. 

Figures 5.1.3 and 5.1.4 show results for 30x30 and 30x60 mm delaminations. 
Comparison of Figs. 5.1.3 and 5.1.4 shows that except for the lowest strains, the 
smaller (circular) delamination actually has a larger Gj. Hence, the growth rate 
based on G j is expected to be larger for the smaller delamination. As pointed out 
earlier, the distributions of Gjj for circular and elliptical delaminations are quite 
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(a) Mode I 



(b) Mode II 


Fig. 5.1.5 Strain-energy release rate distributions for a circular delamination, 
(diameter = 60mm) 



different. Even the location of maximum Gjj is different. The peak values of Gjj 
for the two delaminations are similar. Hence, based on Gjj the growth rate should 
be about the same for both delaminations, but the direction of growth would be 
different. 

5.2 Interpretation of Deformation and Load Transfer 

Section 5.1 presented some numerical results for strain-energy release rates. 
Several trends were observed and discussed. The purpose of this section is to 
present an intuitive explanation of the process by which load is redistributed 
and secondary loads created which lead to instability-related delamination growth 
(IRDG). 

The mechanics of IRDG for the through-width delamination have been de- 
scribed previously in refs. 3, 4, 5. The mechanics of IRDG for the embedded 
delamination will be derived here. Both axisymmetric and uniaxial loading of the 
embedded delamination will be considered (even though axisymmetric loading is 
not considered elsewhere in this thesis). The approach will be to present first G/ 
and Gjj results which illustrate the different behaviors of several configurations 
which exhibit IRDG. Then mechanics arguments will be offered to explain the 
different behaviors. 

5.2.1 Comparison of Behaviors for the Through- Width and 
the Embedded Delamination 

Four cases were analyzed: the through-width delamination, the axisymmet- 
rically loaded embedded circular delamination, the uniaxiaily loaded embedded 
circular delamination, and the uniaxiaily loaded embedded elliptical delamina- 
tion. In all cases the thicknesses H and h were 4 and .4 mm, respectively. The 
through-width delamination was 30mm long. The circular delaminations were 
30mm in diameter. The elliptical delamination was elongated along the y-axis 
and had dimensions of 30x60mm. The same strain range was used for all the 
cases. 
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Figures 5. 2. 1.1 and 5.2. 1.2 show the effect of strain level on Gj and Gjj for 
the through-width delamination and the embedded delamination. The mode III 
component is not shown because it was essentially zero for all the cases considered. 
In general, one would not necessarily expect Gjjj to be zero. Both uniaxial 
and axisymmetric loads were considered for the embedded delamination. For 
the uniaxially loaded embedded delamination, both a circular and an elliptical 
delamination were analyzed. Figs. 5.2.1.1c, 5.2.1. Id, 5.2.1.2c, and 5. 2.1. 2d show 
the variation with strain level at two points along the delamination front of the 
uniaxially loaded embedded delamination: at 0 = 0° and 0 = —90°. See Fig. 2.8.2 
for the definition of 0. 

The variation of Gj with strain level is dramatically different for the four 
cases. For the through-width delamination Gj increases very rapidly after buckling 
occurs (Fig. 5.2.1.1a). After reaching a peak, Gj decreases. For the axisymmetric 
case, Gj increases monotonically. Also, the magnitude of Gj in Fig. 5.2.1.1b is 
much larger than in Fig. 5.2.1.1a. 

Fig. 5.2.1.1c shows the variation of G/ at 0 = -90° for a uniaxially loaded 
embedded delamination. Only the results for the elliptical delamination are 
shown here, since Gj was zero at 0 = —90° for all strain levels for the circular 
delamination. The shape of the curve for the elliptical delamination is similar 
to that for the through-width delamination (Fig. 5.2.1.1a), but the magnitude is 
less. At 6 = 0°, Gj is shown for both the circular and elliptical delaminations 
(Fig. 5. 2.1. Id). Note that the magnitude of Gj is much larger than in Figs. 
5. 2.1. la-5.2. 1.1c. Also, Gj increases rapidly and monotonically with applied strain. 

Fig. 5. 2. 1.2 shows the Gjj variations with strain level. In all cases Gjj 
increased monotonically with strain. Gjj is of the same order of magnitude for 
all the cases except for the uniaxially loaded circular delamination at 0 = —90° 
(Fig. 5.2.1.2c). This contrasts with the very wide range of magnitudes in 
Fig. 5.2.1. 1 for Gj. 
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(c) (d) 


Fig. 5. 2. 1.1 Effect of strain level on Gj. (a) Through-width delamination, 
(b) Embedded delamination with axisymmetric load, (c) Embedded delamination 
with load at B = -90, and (d) Embedded delamination with uniaxial load at d = 0. 
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(c) (d) 


Fig. 5.2. 1.2 Effect of strain level on G//. (a) Through-width delamination, 

(b) Embedded delamination with axisymmetric load, (c) Embedded delamination 
with uniaxial load at 6 - -90, and (d) Embedded delamination with uniaxial load 
at 6 = 0. 
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Even though the strain-energy release rates in all the cases illustrated in 
Figs. 5. 2. 1.1 and 5. 2. 1.2 are a result of local buckling, the variety of behavior 
suggests that there must be variations in the mechanism by which local buckling 
causes strain- energy release rates. The next section will attempt to explain these 
mechanisms. 

5.2.2 Mechanics of IRDG for the Embedded Delamination 

In highly simplified anthropomorphic terms, a strip of the buckled region which 
is parallel to the load direction (strip A in Fig. 5.2.2. l) wants to buckle outward. A 
strip of the buckled region which is perpendicular to the load direction (strip B in 
Fig. 5.2.2. 1) has no desire to deform outward; it is pushed out by strip A. Strip A is 
analogous to the through-width case. The constraint provided by strip B reduces 
Gj for strip A. Conversely, strip A causes high Gj at the ends of strip B when 
strip A pushes strip B outward. Of course, the buckled region is not comprised 
of strips, but this simplified interpretation helps explain the behavior observed. 
The following paragraphs present a more rigorous and detailed discussion of the 
mechanics of instability-related delamination growth (IRDG). 

The through-width delamination will be discussed first. After describing 
the mechanics for the through-width delamination, its close relationship with 
the embedded delamination with axisymmetyric loads will be discussed. Next, 
tractions will be applied to the through-width delamination configuration which 
transform it into a uniaxially loaded embedded delamination. The required 
tractions should give some feel for why the behaviors differ for the embedded 
delamination and the through-width delamination under uniaxial loads. 

The discussion of the through-width delamination can be expedited by first 
transforming this geometrically nonlinear problem into a linear one with nonlin- 
early related loads (ref. 4). Fig. 5.2.2.2a shows a schematic of a laminate with 
a postbuckled through-width delamination. In Fig. 5.2.2.2b the buckled region is 
replaced by the loads Pq and M, the axial load and moment,respectively, in the 
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Fig. 5.2.2. 1 Strips of sublaminate aligned perpendicular and parallel to the load 
direction. 
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Fig. 5.2 .2.2 Nonlinear configuration (a) transformed into linear configuration (d) 
with two nonlinearly related loads (Pq — Pp) and M. 


buckled region where it is cut. The total applied load Pj is equal to P B + P c . The 
load system in Fig. 5.2.2.2c, (which is the same as Fig. 5.2.2.2b) can be divided 
into the two load systems shown in Figs. 5. 2. 2. 2d and 5.2.2.2e. Because Pq and 
P B are calculated using rule of mixtures, the load system in Fig. 5.2.2.2e causes 
a uniform axial strain state and no interlaminar stresses. Accordingly, for strain- 
energy release rate calculation, the configuration in Fig. 5.2.2.2d is equivalent to 
the original configuration (Fig.5.2.2.2a). The moment M opens the delamination, 
contributing to G/. It also contributes to some G// (ref. 4). Also, the load 
(p c _ p D ) contributes to G/j. In addition, because of the offset of the line of 
action of Pc ~ Pd relative to the delamination, this force creates a moment which 
tends to close the delamination and reduce the Gj component caused by M. The 
result of the competing mechanisms are strain-energy release rate variations like 
that in Fig. 5a and 5.2.1.2a. The mode I strain-energy release rate first increases 
very rapidly with increasing strain and then decreases to zero. The mode II strain- 
energy release rate increases monotonically with applied load, since both M and 
(Pc- Pd) contribute to G//. 

The axisymmetrically loaded circular delamination is very similar to the 
through-width delamination. In fact, the schematics in Fig. 5.2.2.2 are applicable 
if the forces are replaced by forces per unit length. The load in a column, Pd, is 
essentially constant after the applied strain is increased beyond the buckling strain, 
but the load in an axisymetrically loaded plate continues to increase significantly 
after buckling. The load Pc increases linearly with the applied load. Hence, (Pc~ 
p D ) is large for the through-width delamination but is relatively small for the 
axisymmetric case. As a result, there is little attentuation of the effects of M by 
(P c - P D ) for the axisymmetric case. Consequently, G/ is much larger for the 
axisymmetric case (Compare Figs. 5.2.1.1a and 5.2.1.1b). 

Now the uniaxially loaded embedded delamination will be considered. Figs. 
5.2.2. 3-5. 2.2.6 illustrate the transformation of a through-width delamination 
(Fig. 5.2.2. 3) into an embedded delamination (Fig. 5.2.2.6). The letters A through 
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Fig. 5. 2. 2. 4 Slice of laminate showing tractions required to perform transformation 
by closing buckled part of boundary BF. 
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Fig. 5. 2.2. 5 Slice of laminate after tractions have closed boundary. 
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G are added to aid the discussion. They are not related to the letters in Fig. 5. 2. 2. 2. 
To expedite the discussion, the embedded delamination will be assumed to be rect- 
angular. Tractions are required to close the buckled part of region AEFB. These 
will only be non-zero near the new delamination front and are, in fact, the inter- 
laminar stresses. 

Fig. 5. 2. 2. 4 shows a slice removed from the laminate. The figure shows the 
forces required to close the buckled part of the boundary BF. These forces are 
generated when the region AEFB is closed. These forces indicate some of the 
interaction of regions AEFB with BFGC. There are in-plane forces F y , transverse 
forces F z , and a moment M y . Fig. 5.2.2.5 shows the same slice after the forces 
have closed the buckled part of BF. 

The moment M y would operate in the direction indicated in Fig. 5. 2. 2.4 based 
on the curvature in Fig. 5. 2. 2. 5. Likewise, the transverse force, F z , would be 
expected to act downward to help close the delamination front. The force, F y , is 
a result of two things : transverse deflection and Poisson’s ratio. When transverse 
deflection occurs, the length of a line from A to D (Fig. 5.2. 2.6) increases, hence 
the buckled region must be stretched in the y-direction. When the laminate is 
compressed in the x- direction, it expands in the y-direction due to Poisson’s effect. 
If the base laminate has a larger Poisson’s ratio than the sublaminate, then a force 
F y is required to enforce compatibility when the buckled part of BF is closed. The 
sign of F y due to Poisson’s ratio would depend on the relative magnitudes of the 
Poisson’s ratios. Differences in Poisson’s ratio were not considered in this study. 

The magnitude of F y should be related to the in-plane stiffness in the y- 
direction. The magnitude of M y should be related to the flexural stiffness in 
the y-direction. The effect of material properties on F z is not as straight-forward, 
so no prediction will be offered. 

The dimensions of the embedded delamination should affect the forces and 
moment. For the same transverse deflection, the curvature in the y-direction is 
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less for a larger dimension b , so M y should decrease with increasing b. Also, the 
strain in the y-direction due to the increased length A-D (Fig. 5. 2. 2. 6) would be 
less for larger b. Hence, for the same transverse deflection F y should decrease with 
increasing b. 

The moment M y should contribute primarily to Gj, but it also contributes to 
Gjj. The force F y should contribute to Gjj and reduce Gj. The reduction is due 
to the offset between the delamination plane and the middle of the sublaminate. 
This offset causes a moment relative to the delamination plane which is opposite 
to M y . Based on the large G/ in Fig. 5. 2.1. Id, the opening effects of M y must 
dominate the closing effects of F y . 

The original configuration with a through-width delamination (Fig. 5. 2. 2. 3) 
had some distribution of G/ and Gjj along x= a. The application of the forces 
F y , F t , and the moment M y changes the load flow significantly. Closing the 
ends of the through-width delamination (area AEFB) contributes a compressive 
component of interlaminar stresses along the front x=a, thus reducing G/. In fact, 
because of this reduction, G/ was zero for all strain levels at 6 = —90° for the 
circular delamination under uniaxial loads. For the elliptical delamination, the 
magnitude of G/ at 6 = -90° was less than for the through-width delamination 
for the same reason (Fig. 5.2.1.1a and 5.2.1.1c). 

Based on the preceeding discussion, one would expect the behavior of a 
through-width delamination and the uniaxially loaded embedded delamination 
to have some similarity, but probably more differences. Fig. 5. 2. 1.1 illustrates this 
very well. For an embedded delamination highly elongated perpendicular to the 
load direction, one would expect the behavior to be like that for the through-width 
delamination. No highly elongated delaminations were examined in this study, but 
even the 1:2 aspect ratio ellipse has a Gj variation at 8 = —90 (Fig. 5. 2.1.1c) 
which is very similar to that for a through-width delamination (Fig. 5.2.1.1a) . 
But this similarity has little importance for this case, since the G/ was very much 
larger at 9 = 0° (Fig. 5. 2.1. Id). 
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5.3 Comparison of Full 3D. Thin-film 3D. and Thin-film Plate Analyses 
For an Embedded Delamination 


Geometrically nonlinear 3D finite element analysis is inherently very expensive. 
Hence, there is considerable motivation to use simplified and less rigorous (and less 
expensive) analysis techniques. An obvious approximation is to assume that the 
base laminate is so much thicker and stiffer than the sublaminate, that the base 
laminate strains are independent of how the sublaminate deforms. For most of the 
configurations analyzed herein, the base laminate deformation can be expressed in 
closed form. Because of the base laminate dominance, only the sublaminate needs 
to be analyzed. This is referred to as a thin-film analysis. Boundary conditions 
consist of displacements prescribed to impose compatability of the sublaminate 
with the known base laminate deformations. This thin-film assumption can be 
used with 3D or with plate analysis. 

The thin-film assumption has been used with plate analysis by several re- 
searchers, as mentioned in Chapter 1. However, only ref. 21 gives distributions of 
total strain-energy release rate. This section will compare results for a square and 
a rectangular delamination from ref. 21 with 3D solutions.The sublaminates are 
assumed to have an initial waviness of a sinusoidal shape, as given by eqn. 5.3.1. 


.05, 7rx w Try . , * 

I "initial = —(1 + cos — )(1 + cos — ) mm (5.3.1) 

The sublaminate and base laminate thicknesses were .51 mm and 5.1 mm, 
respectively. The entire laminate is assumed to be isotropic with a Young’s 
modulus of 53.3 GPa and a Poisson’s ratio of .31. Both full 3D and thin-film 
3D results will be presented. 

Fig. 5.3.1 shows the distribution of Gj along the delamination front for the 
square delamination (25.4 x 25.4 mm). Three strain levels were considered. Results 
are shown for full 3D, thin-film 3D, and thin-film plate analysis. There is good 
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Fig. 5.3.1 Total strain-energy release rate distribution for a square delamination 
calculated using 3D, thin-film 3D, and plate analysis. (2a x 2b = 25.4 x 25.4mm) 
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agreement between the three analyses for all three strain levels. Fig. 5.3.2 shows 
the distribution of Gj along the delamination front for a rectangular delamination 
(25.4 x 50.8 mm). The agreement between the three analyses is not as good as for 
the square delamination, but it still is fairly good. 

For the cases compared, thin-film plate analysis gave reasonably good predic- 
tions of Gj. When one considers the high sensitivity of Gj 1 to parameters such 
as sublaminate thickness, delamination size, and initial imperfection and the un- 
certainty of knowing these parameters precisely, the differences between the plate 
analysis and 3D analysis appear negligible. Of greater importance is the limitation 
that plate analysis can only calculate Gj, not the components. A hybrid analysis 
which uses 3D analysis near the crack front and plate analysis elsewhere appears 
to be a good alternative. Such a technique was used in ref. 4 for the through-width 
delamination. 

The results in Figs. 5.3.1 and 5.3.2 also serve a secondary purpose. Totally 
different analyses were used to obtain the results in ref. 20 and the 3D results. 
Also, the method for calculating strain- energy release rates were different. But the 
agreement is good. Hence, the results give additional validation to both the virtual 
crack closure technique presented in ref. 20 and to the 3D analysis developed 
herein. 

5.4 Effect of Contact Constraints 

Earlier it was shown that for certain combinations of delamination size and 
strain level, closure occurs over part of the delamination front. Unless constraints 
are imposed to prevent interpenetration, overlap of the crack faces will occur in 
an analysis. The results presented in section 5.1 are based on an analysis which 
allowed interpenetration of the crack faces. This section examines the effects of 
including contact constraints on the calculated distribution of strain-energy release 
rates. 

The configuration analyzed is a “homogeneous” quasi-isotropic laminate with 
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Fig. 5.3.2 Total strain-energy release rate distribution for a rectangular de- 
lamination calculated using 3D, thin-film 3D, and plate analysis. (2a x 26 = 
25.4 x 50.8mm) ( e x = —.0043) 
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a circular delamination (2a x 2b = 30 x 30 mm). The thickness of the sublaminate 
is .4 mm and the base laminate is 4.0 mm. For the circular delamination, closure 
occurred near $ = —90° for all of the strain levels considered. The amount of 
closure at 90° relative to the amount of opening elsewhere along the front increased 
with the ratio of applied strain to the bifurcation buckling strain. Hence, only large 
strains will be considered here: -.005 and -.02 . These are about 3 to 11 times the 
bifurcation buckling strain. Of course, -.02 strain is unrealistically large for current 
composites. But larger delaminations or thinner sublaminates can be subjected 
to strains more than 11 times the buckling strain. The -.02 strain case should 
simply be considered a case where the applied strain is quite large compared to 
the buckling strain. 

Fig. 5.4.1 and 5.4.2 show the deformation of the sublaminate for the two 
strains before and after imposing contact constraints. Only the top surface of 
the sublaminate is plotted. The displacements were scaled up by a factor of 10 
and 2.5 for the cases e x = —.005 and -.02, respectively. The interpenetration is 
much more severe for the larger strain. 

Fig. 5.4.3 shows Gj and Gjj for a strain of -.005. For comparison, results 
with and without contact constraints are shown. Because of the approximations 
mentioned earlier in imposing the constraints, Gj is not computed to be identically 
zero in the contact region. But it is now negligibly small. In the region where the 
delamination front is open, Gj is not very sensitive to whether or not contact 
constraints are imposed. In absolute terms, Gjj is also not affected much. 
Fig. 5.4.4 shows Gj and Gjj for a strain of -.02. Since the overlap (before 
constraints are imposed) is much greater for the larger strain, one would expect 
a larger effect on Gj and Gjj. Comparison of Figs. 5.4.3 and 5.4.4 shows this is 
the case. But the trends are the same for both strain levels. That is, imposition 
of the constraints reduces Gj to a negligible magnitude in the contact region and 
increases Gj in the non-contact region. Gjj is increased along the entire front. 

For a strain level of -.005, little error is incurred by not imposing contact 
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(a) without contact constraints 


Fig. 5.4.1 Deformation of top surface of postbuckled sublaminate with and without 
contact constraints. Displacements multiplied by 10. (er x = -.005) 
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(b) with contact constraints 


Fig. 5.4.1, Concluded. 
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(a) without contact constraints 


Fig. 5.4.2 Deformation of top surface of postbuckled sublaminate with and without 
contact constraints. Displacements multiplied by 2.5. (e* = — -02) 
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(b) with contact constraints 


Fig. 5.4.2, Concluded. 
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Fig. 5.4.3 Effect of contact constraints on Gj and Gjj distributions (e x = —.005). 
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Fig. 5.4.4 Effect of contact constraints on Gj and G// distributions (e x — 


-. 02 ). 
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constraints and just recognizing that G/ is zero in the contact region. For a 
strain of -.02, significant errors occur when contact constraints are not imposed. 
Since there is a significant computational cost associated with including contact 
constraints, constraints probably should not be imposed unless the overlap is quite 
extensive. 

5.5 Strain-Energy Release Rates For an Edge Delamination 

This section will discuss the results of a limited parametric study of homoge- 
neous quasi-isotropic laminates containing a postbuckled edge delamination. The 
parameters varied were strain level, delamination shape, and delamination size. 
The effects of these parameters on deformation and strain-energy release rate were 
determined. 

Figure 5.5.1 shows the variation of G/, G//, and G<p along the delamination 
front for a semi-circular edge delamination of radius 15mm. The strain-energy 
release rates increase rapidly with increased strain. There is a large variation of 
the strain-energy release rates along the delamination front. The maximum value 
of G/ is a little larger than the maximum value of G// for each strain level. Note 
that there is no overlap region for this configuration, in contrast to the response 
for the circular embedded delamination (see Fig. 5.1.3). The maximum G/ and 
Gn occur at different locations for the edge delamination, which is also different 
than the circular embedded delamination behavior. In fact, the semi-circular 
edge delamination behaves somewhat like an elliptical embedded delamination. 
Fig. 5.1.4 shows that an elliptical embedded delamination exhibits little overlap 
and has maximum values of G/ and G// at the same locations that the edge 
delamination has its maximum values. However, G// increased monitonically 
with perimeter coordinate S for the elliptical embedded delamination. For the 
edge delamination G// decreases and then increases with S. For the elliptical 
embedded delamination G'p was almost constant for some strain levels. The Gf 
is far from constant for any strain level for the edge delamination. The predicted 
direction of delamination growth depends on the assumed growth criterion, since 
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the peak G/ and G// occur at different locations. However, growth should occur 
preferentially either near the x- or y-axes, and not in between. 

Figure 5.5.2 shows the variation of G/, G//, and Gj along the delamination 
front for a semi-elliptical edge delamination. The maximum values of Gj and G// 
still occur at different locations. Now the maximum value of G// is much larger 
than that for G/. The predicted direction of growth is still dependent on the 
assumed growth criterion. However, the semi-elliptical is more likely to grow in 
the load direction than is the semi-circular delamination. 


Figure 5.5.3 shows the variation of G/, G//, and Gy along the delamination 
front for a semi-circular edge delamination of radius 30mm. The distributions are 
similar to the distributions in Fig. 5.5.1 for the smaller semi-circular delamination. 
The magnitudes of the strain-energy release rates are larger for the larger semi- 
circular delamination. Note that in Fig. 5.5.3 there is a small amount of overlap, 
which did not occur for the smaller delamination. This is because the ratio of the 
applied strain to the buckling strain is larger for the larger delamination. 
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CHAPTER 6 

DESCRIPTION of SPECIMENS and EXPERIMENTAL PROCEDURES 

This chapter describes the various specimens which were tested and the 
procedures used to conduct the tests. There are two subsections in this chapter . 
1. Specimen configurations and 2. Measurement of deformation and damage. 

6.l Specimen Configurations 

There were three basic specimen configurations: the transversely loaded 

plate, a laminate with an embedded delamination, and a laminate with an edge 
delamination. 

Fig. 6.1.1 shows a schematic of the transversely loaded plate configuration. A 
thin sheet of steel was bonded to a thick aluminum plate with a room temperature 
cure adhesive. The aluminum plate had a 50.5mm diameter through-hole. A single 
transverse load was applied to the center of the steel using a 12.7mm diameter 
indenter. Fig. 6.1.2 shows the specimen and the test fixture. This set-up is the 
same as that used in ref. 44. Direct-current differential transformers were used 
to monitor load point deflections. The purpose of these tests were to obtain load 
versus deflection. No debond measurements were performed on these specimens. 

Several types of specimens with embedded delaminations were fabricated. One 
of the specimens consisted of spring steel bonded to aluminum. Fig. 6.1.3 shows a 
schematic of the specimen. This specimen type was tested to provide postbuckling 
deformation data for a laminate without stacking sequence effects. As indicated 
in the figure, EA934NA room-temperature cure adhesive was used for bonding. 
A teflon insert was used to prevent bonding over a circular region, resulting in a 
simulated delamination. 

Another type of specimen with an embedded delamination consisted of a thin 
AS4/PEEK laminate with a delamination. The stacking sequence was either 
(0/90/90/0), or (90/0/0/90),. Fig. 6.1.4 shows a photograph of a typical laminate. 
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Fig. 6.1.1 Schematic of transversely loaded plate specimen. 
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Fig. 6.1.3 Steel/ A1 specimen with embedded delamination. 




As shown in the figure, the laminates were badly warped. The PEEK laminate was 
bonded to an aluminum sheet to increase resistance to global buckling. Fig. 6.1.5 
shows a schematic of the specimen. A double layer of thin Kapton film (thickness 
of each layer= ,013mm) was used to provide an initial delamination of 30, 40 or 
60mm diameter between the fourth and fifth plies. The double layer of Kapton 
was used because Kapton might bond with the PEEK, but it should not bond to 
another piece of Kapton. Apparently the high temperature and pressure used in 
processing the PEEK laminates caused a small amount of bonding. Even a small 
amount of bonding was sufficient to prevent local buckling of the delamination 
region prior to global collapse. Hence, a small block of aluminum was bonded to 
the sublaminate using a low strength bond, as illustrated in Fig. 6.1.6, and small 
forces were applied by hand to loosen the Kapton bond. Because of this poorly 
quantified “preconditioning” and the bad warpage problem, the specimens were 
only used for qualitative measurements of delamination growth. 

Twenty-four ply composite specimens with embedded delaminations were 
fabricated from IM7/8551-7 prepreg. The stacking sequences were (0/90/90/0) 6 
and (90/0/0/90)6- A double layer of Kapton film was used to provide an initial 
delamination between the fourth and fifth plies. Delamination sizes were 30, 40, 
or 60mm. Unlike the PEEK laminates, the Kapton film did not exhibit spurious 
sticking in these specimens. 

Several of the PEEK specimens described above were sliced longitudinally to 
obtain laminates with semi-circular edge delaminations. 

The specimens with embedded or edge delaminations were tested in compres- 
sion. Steel guide plates were used to prevent global buckling. Most of the tests 
were conducted using a solid guide plate on one side of the specimen and a plate 
with a 82mm window on the other side. A few tests were performed using plates 
with windows on both sides. The guide plates are shown in Fig. 6.1.7. The guide 
plates were 19mm thick. For specimens with an embedded delamination, the 
window was centered over the delaminated region. For specimens with an edge 
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Fig. 6.1.6 Breaking Kapton to Kapton bond with glued on block. 
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delamination, the guide plates were offset to one side to provide extra support. 

6.2 Measurement of Damage 

Several techniques were used for detecting damage in the compression spec- 
imens. Because delamination growth resulted in a larger buckled region, visual 
inspection during the loading proved to be fairly accurate for detecting when de- 
lamination growth initiated. A more refined technique was measurement of the 
deformed shape of the buckled region. The delamination boundary corresponded 
to the location where transverse displacements were negligible. The measurements 
were obtained using the fixture illustrated in Fig. 6.2.1. This fixture is very similar 
to that described in ref. 45. Basically, the fixture consists of a machinist’s vise 
with three DCDT’s. One of the DCDT’s measures the transverse deflection of the 
buckled region and the other two indicate where the transverse deflection is being 
measured. In some cases, two DCDT’s were used to monitor transverse displace- 
ments; one DCDT monitored back face displacement and the other monitored front 
face (i.e. where the buckling occurred) displacement. By subtracting the outputs 
of the two transducers, a change in thickness could be measured. Essentially all 
change in thickness would be due to postbuckling displacements. This technique 
using two DCDT’s was less sensitive to artifacts due to slight misalignment of the 
specimen and the machinist’s vise. 

Another technique for measuring damage was applying an X-ray opaque 
dye penetrant and then taking an X-ray. The dye highlighted ply cracks and 
delaminations. However, if little or no surface damage occurred, the dye penetrant 
was not able to reach the internal damage area. In those cases, after all testing 
of those specimens were completed, a small hole was drilled in the center of the 
specimen. Then, the dye could be applied internally. 
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Fig. 6.2.1 Deflection measurement fixture. 


Some specimens were sectioned using a diamond saw and polished. The cross 
sections were examined using light microscopy. This technique was particularly 
useful for determining the distribution of cracking through the thickness and 
details of the types of damage. 
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CHAPTER 7 

COMPARISON of ANALYTICAL and EXPERIMENTAL RESULTS 


The primary purpose of this chapter is to describe the results of a combined 
experimental and analytical study of instability-related delamination growth. 
Specimens were designed based on the trends observed in the analytical study 
discussed earlier. Based on the analytical study differences in the buckling strain, 
the shape of the postbuckled region, the strain at which delamination growth would 
occur, and the direction of growth were expected. The experimental program was 
conducted to determine whether the expectations based on analysis were correct. 
Also, it was hoped that a strain-energy release rate parameter would be identified 
which could be used to quantitatively predict the onset of delamination growth. 

Two composite material systems were considered: AS4-PEEK and IM7/8551-7. 
Because of the processing difficulties described earlier in Chapter 6, only quali- 
tative comparisons of analysis and experiments will be presented for the PEEK 
specimens. 

A secondary purpose of this chapter is to present results for two “check cases”: 
the transversely loaded plate and the steel bonded to aluminum compression 
specimen. These configurations were examined because they were expected to 
behave in a fairly predictable fashion. Hence, they provided a good starting point 
for comparison of analytical and experimental results. These check cases will be 
presented first, £nd then the primary study will be discussed. 

7.1 Check Cases 

Two types of experiments were conducted just to help verify the finite element 
program NONLIN3D. One of the experiments involved central transverse loading 
of a circular plate (see Fig. 6.1.2). The other test consisted of spring steel bonded to 
aluminum and loaded in compression to cause local postbuckling of a delaminated 
region (see Fig. 6.1.3). The results for these two check cases are described in this 
section. 
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Figure 7.1.1 shows experimentally measured load point deflection versus ap- 
plied load for the transversely loaded plate. Results are shown for two specimens. 
Each specimen was loaded and unloaded more than once, so there are multiple 
data for each load level. Also shown are the results from NONLIN3D. The fig- 
ure shows that NONLIN3D predicts the deflection quite well, even though the 
response is highly nonlinear. 

Figure 7.1.2 shows results from the steel/ A1 postbuckling tests. This figure 
shows the deformation of the specimens and the finite element models along the 
planes x = 0 and y = 0. The finite element analysis was performed at a slightly 
different strain level than was present in the test. Ideally, the same strain would 
have been used for both the analysis and the experiments. However, the deflection 
increases rapidly with applied strain and even small imperfections can shift the 
load versus deflection curve. To expedite the comparison of the shape of the 
buckled region with the prediction, the strain for the finite element analysis was 
chosen so that the peak deflection would be about the same as was observed in 
the experiment. 

Figure 7.1.2 shows the deformation for a specimen with a 60mm diameter 
delamination. The predicted and observed deformed shape agree quite well. Note 
that the deformation is much different along the x = 0 and y = 0 planes. Along 
the y = 0 plane the large deformation is restricted to a much narrower region than 
along the x = 0 plane. This restricted deformation caused the delaminated region 
to appear elliptical during the test, even though the delamination was known to 
be circular. 

7.2 Qualitative Study of AS4/PEEK Laminates 

As described earlier in the experimental procedures section, there were prob- 
lems with the fabrication of the AS4/PEEK specimens. Consequently, only a qual- 
itative study was performed on this material system. Observations were made of 
the shape of the buckled region, the location of delamination growth, the direction 
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Fig. 7.1.1 Comparison of measured and predicted load point displacement for a 
transversely loaded circular plate. 
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(a) x = 0 plane 


Fig. 7.1.2 Deformation of the x = 0 and y = 0 planes in a steel/AL specimen with 
a postbuckled 60mm diameter delamination. 
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Experiment 



Fig. 7.1.2, Concluded 


of growth, and ply cracking. Fatigue loads were used to cause delamination growth 
in most of the tests because of the high resistance of PEEK to static delamination 
growth. A few specimens were tested statically. However, effects due to problems 
with specimen fabrication prevented isolation of fatigue versus static behavior. 
Laminates with either an embedded or an edge delamination were tested. Two 
sublaminate stacking sequences for the sublaminate were considered: (0/90/90/0) 
and (90/0/0/90). The embedded delamination will be discussed first, then the 
edge delamination. 

Fig. 7.2.1 shows sketches of two specimens (which had a 60mm diameter initial 
delamination). The boundary of the initial delamination is indicated by dashed 
lines. The region which actually buckled outward is indicated by solid lines. This 
region was determined by visual inspection. Only part of the delaminated region 
buckles outward for both sublaminates. Also, the buckled region is much narrower 
for the (90/0/0/90) sublaminate than for the (0/90/90/0) sublaminate. The same 
trends were observed for 40mm delaminations. 

Finite element analysis was performed before the tests were conducted. In 
fact, the stacking sequences for the experiments were selected because the analysis 
predicted that the behaviors of the two laminates should be different. These pre- 
liminary analyses were performed for a 30mm delamination. Since problems with 
fabrication precluded quantitative comparisons, the analyses were not repeated for 
the actual delamination sizes. Fig. 7.2.2 shows the y=0 plane for deformed finite 
element models for the two stacking sequences. These particular models were for 
a 30mm delamination, but the trends are independent of delamination size. Note 
that the analysis prediction agrees with the experiments, i.e. that the buckled 
region should be narrower for the 90/0/0/90 sublaminate than for the 0/90/90/0 
sublaminate. 

Fig. 7.2.3 shows strain-energy release rate distributions for a 30mm delami- 
nation for the two stacking sequences. Both Gj and Gjj are maximum at S=0 
and are much less elsewhere. These results suggest that delamination growth 
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Fig. 7.2.1 Outline of postbuckled region for 60mm delamination in two AS4/PEEK 
laminates. 


w 0 = .73 mm 



(a) (0/90/90/0) sublaminate 


w 0 = .70 mm 



(b) (90/0/0/90) sublaminate 


Fig. 7.2.2 Deformation of the y = 0 plane for two AS4/PEEK laminates. Diameter 
of delamination is 30mm. 
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(a) (0/90/90/0) 9ublaminate 

Fig. 7.2.3 Strain-energy release rate distributions for two AS4/PEEK sublami- 
nates. Delamination diameter is 30mm. 
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should initiate in a small area and should propagate transverse to the load di- 
rection. Figure 7.2.4 shows X-rays for the two stacking sequences. The observed 
initial delamination growth is very localized and the growth is transverse to the 
load direction. Of course, after a small amount of growth the analytical results 
in Fig. 7.2.3 are not applicable, since the configuration has changed considerably. 
The X-rays show matrix cracking, which is apparently caused by the flexure during 
postbuckling. The delamination growth was usually very rapid. Sometimes there 
was essentially instantaneous propagation to the guide plates, but not always. 

Stacking sequence effects were quite obvious for the edge delammation speci- 
mens also. It was difficult to obtain delamination growth for the 30mm diameter 
delamination for the (0/90/90/0) sublaminate. Except for one case, global failure 
occurred simultaneously with delamination growth. In the one case that growth 
could be observed before global failure, only a small amount of growth occurred 
and approximately the same amount of growth occurred in both the load direction 
and transverse to the load direction. Fig. 7.2.5 shows an X-ray of a laminate with 
a 20mm radius initial delamination and a (0/90/90/0) sublaminate. In this case, 
the dominant delamination growth is definitely in the load direction. Also, when- 
ever catastrophic growth occurred, there was significant delamination growth in 

the load direction. 

Only 30mm delaminations were examined for the (90/0/0/90) sublammmate. 
The X-ray in Fig. 7.2.6 shows that for this case the delamination growth is 
definitely transverse to the load direction. Even when extensive delamination 
growth occurred, it was transverse to the load direction. Matrix cracking was 
more obvious for the 90/0/0/90 sublaminate (Fig. 7.2.6) than for the 0/90/90/0 
sublaminate (Fig. 7.2.5). 

Figures 7.2.7 and 7.2.8 show strain-energy release rate distributions for two 
specimens with edge delaminations. The analytical results are for approximately 
the same range of strains that were used in the tests. Results are shown for three 

strain levels. 


147 





148 


ORIGINAL PAGE IS 
OF POOR QUALITY 


Fig. 7.2.4 Radiographs of two peek laminates with 60mm diameter initial 
dclarnination. 
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(b) (90/0/0/90) Sublaminate 
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Fig. 7.2.5 Radiograph of PEEK laminate with an initial semi-circular edge 
delamination of radius 20mm. The sublaminate is (0/90/90/0). 
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Fig. 7.2.7 Strain-energy release rate distributions for (0/90/90/0) AS4/PEEK 
sublaminate with a semi-circular edge delamination. 
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Fig. 7.2.8 Strain-energy release rate distribution for a (90/0/0/90) AS4/PEEK 
sublaminate with a semi-circular edge delamination. 




The strain-energy release rate distributions for a 40mm edge delamination and 
(0/90/90/0) sublaminate are shown in figure 7.2.7. Based on G/ one would expect 
growth transverse to the load direction. But based on G//, one would expect 
growth along the load direction. Without a mixed-mode delamination growth 
criterion, it is not possible to predict the direction of delamination growth. The 
experiments usually exhibited dominant growth in the load direction. The figure 
shows that the peak G// is about twice the peak Gj. Hence, it appears that the 
large magnitude of G// was sufficient to cause dominant growth to occur in the 
load direction. 

The strain-energy release rate distributions for a 30mm delamination and a 
(90/0/0/90) sublaminate are shown in Fig. 7.2.8. The peak G// component is 
much smaller than the peak G/, and would not be expected to play much of a role. 
Based on G/ one would clearly expect dominant delamination growth to occur 
transverse to the load direction. This agrees with the experimental observations. 

In summary, there are significant qualitative agreements between the observed 
and predicted behaviors. 

7.3 Quantitative Study of IM7/8551-7 Laminates 

A variety of measurements and observations were made for IM7/8551-7 lami- 
nates with an embedded delamination. As was the case for PEEK, two stacking 
sequences were used: (0/90/90/0)6 and (90/0/0/90)6 • The initial delamination 
was located between the fourth and fifth plies. No edge delamination tests were 
performed. Because there were no major manufacturing induced artifacts (like 
the Kapton sticking problem in the PEEK laminates), quantitative comparisons 
were made between analysis and experiments for postbuckling deformations and 
initiation of delamination growth. Also, the stability of delamination growth was 
monitored. A combination of techniques were used to determine the extent and 
types of damage. These techniques included measurement of deflections, X-ray, 
ultrasonic C-scans, and sectioning followed by light microscopy. 
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Figures 7.3. 1-7. 3.3 show load versus deflection for the six laminate configura- 
tions considered (2 stacking sequences and three delamination sizes). The deflec- 
tion is that in the center of the delaminated region. For each delamination size 
the agreement between the observed and predicted load versus deflection is fair for 
the (0/90/90/0) stacking sequence. For the (90/0/0/90) stacking sequence, the 
agreement is fair for the 30 mm delamination. For the 40 and 60mm delaminations, 
there is a large discrepancy, except at the lower strain levels. This divergence at 
the larger strains may be due in part to the overlapping of the delamination faces, 
which occurs in the analysis, but not in reality. The significance of the overlap 
tends to increase as the ratio of the applied strain to the bifurcation buckling 
strain increases. This ratio increases with strain for a fixed delamination size and 
with delamination size for a fixed strain. This may explain the better agreement 
at lower strains and for smaller delamination sizes. 

Figures 7.3.4 and 7.3.5 show the deformed shape of the y=0 plane for the 
two stacking sequences and a 60mm delamination as measured during tests and 
as predicted from analysis. The strain for the tests was -.0026. The analytical 
results are for a strain of -.003. This small difference in strain is not important 
for the comparisons of deformed shape. There is considerable difference in the 
deformation for the two laminates for the y=0 plane. The buckled area is 
noticeably wider for the (0/90/90/0) laminate than it is for the (90/0/0/90) 
laminate. The deformation for the x=0 plane is not much different for the two 
laminates. The analysis predicted interpenetration of the delamination faces (since 
contact constraints were not imposed). The measured and predicted deformed 
shape agree reasonably well except, of course, where the analysis predicted 
interpenetration. 

Strain-energy release rate distributions were calculated for the six configu- 
rations for a range of strains which bracketed the strain at which delamination 
growth from the Kapton insert occurred. In all cases the G/// component was 
negligible. Figures 7.3.6-7.3.11 show the results of the analysis. There is one figure 
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Fig. 7.3.1 Measured and predicted deflection versus e z for two IM7/8551-7 
laminates. Diameter of initial delamination is 30mm. 
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Fig. 7.3.2 Measured and predicted deflection versus e x for two IM7/8551 
laminates. Diameter of initial delamination is 40mm. 






Fig. 7.3.3 Measured and predicted deflection versus e z for two IM7/8551-7 
laminates. Diameter of initial delamination is 60mm. 
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Fig. 7.3.4 Deformation of the y=0 plane for a (0/90/90/0) sublaminate. Initial 
delamination diameter was 60mm. 
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Fig. 7.3.5 Deformation of the y=0 plane for a (90/0/0/90) sublaminate. Initial 
delamination diameter was 60mm. 



Curve Strain 




Fig. 7.3.6 Strain-energy release rate distributions and maximum strain-energy 
release rate versus s*. The material is IM7/8551-7 and the sublaminate is 30mm 
in diameter and has a stacking sequence of (90/0/0/90). 
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Curve Strain 




Fig. 7.3.7 Strain-energy release rate distributions and maximum strain-energy 
release rate versus fx- The material is IM7/8551-7 and the sublaminate is 30mm 
in diameter and has a stacking sequence of (0/90/90/0). 162 
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Fig. 7.3.8 Strain-energy release rate distributions and maximum strain-energy 
release rate versus s x - The material is IM7/8551-7 and the sublaminate is 40mm 
in diameter and has a stacking sequence of (90/0/0/90). 




Curve Strain 




Fig. 7.3.9 Strain-energy release rate distributions and maximum strain-energy 
release rate versus The material is IM7/8551-7 and the sublaminate is 40mm 
in diameter and has a stacking sequence of (0/90/90/0). 


164 



Curve Strain 




Fig. 7.3.10 Strain-energy release rate distributions and maximum strain-energy 
release rate versus e*. The material is IM7/8551-7 and the sublaminate is 60mm 
in diameter and has a stacking sequence of (00/0/0/90). 
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Fig. 7.3.11 Strain-energy release rate distributions and maximum strain-energy 
release rate versus e x . The material is IM7/8551-7 and the sublaminate is 60mm 
in diameter and has a stacking sequence of (0/90/90/0). 




for each configuration. Each figure shows the variation of G/, G//, and Gj along 
the delamination front. Also shown is the variation of the peak magnitude of G/, 
G//, and Gj 1 with applied strain level. In all cases there is a large variation of the 
strain-energy release rates along the delamination front. The peak magnitudes oc- 
curred at 5 = 0 and the magnitudes were small except close to 5 = 0. Obviously, 
one would expect delamination growth to occur preferentially near 5=0, which 
corresponds to growth transverse to the load direction. During the tests, one could 
visually determine that the delamination growth was localized near 5=0. Fig- 
ures 7.3.12-7.3.15 show X-rays of several specimens after appreciable delamination 
growth had occurred. In these particular specimens there was sufficient through- 
the-thickness damage to permit the dye penetrant to reach the interior. These 
figures show that the analysis was correct in predicting preferential growth near 
S=0. The X-rays show that significant ply cracking accompanied the growth. Ply 
cracking occurred both in 0° and 90° plies. The large horizontal cracks in Figs. 
7.3.12 and 7.2.14 could be seen on the specimen surface with the naked eye. These 
cracks formed after delamination growth had occurred. The tests were stopped 
when the large horizontal cracks formed. 

Several specimens were sectioned and polished after testing. Micrographs of 
these specimens are presented in Figs. 7.3.16-7.3.18. Fig. 7.3.16 shows two 
cross-sectional views of a 0/90/90/0 sublaminate which had an initial diameter 
of 40mm. During the test of this specimen, the delamination growth occurred 
followed (after a" few minutes at constant load) by the formation of a horizontal 
crease in the sublaminate. There was no visible surface damage. Fig. 7.3.16a 
shows that the crease was a result of fiber microbuckling of the interior surface ply 
of the sublaminate. Fig. 7.3.16b shows a micrograph taken near the delamination 
front. Recall that the Kapton implant was located between the fourth and fifth 
plies. The delamination in Fig. 7. 3. 16b is between the third and fourth plies, the 
dpl a.m ina.tion has switched interfaces by way of a matrix crack in the fourth ply. 
This figure also shows angled matrix cracks in the 90° plies. The angled orientation 
was typical in all of the specimens examined. 
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delamination of diameter 30mm. The sublaminate is (90/0/0/90). 






Fig. 7.3.14 Radiograph of IM7/8551-7 laminate with an initial circular embedded 
delamination of diameter 60mm. The sublaminate is (0/90/90/0). 
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(a) View l 


Fig. 7.3.17 Cross-section of IM7/8551-7 laminate with an (90/0/0/90) sublaminate 
and an initial delamination diameter of 60mm. 
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Fig. 7.3.18 Cross-section of IM7/8551-7 laminate with a (90/0/0/90) sublaminate 
and an initial delamination diameter of 40mm. 











Fig. 7.3.17 shows two cross-sections of a 90/0/0/90 sublaminate which had an 
initial di am eter of 60mm. Fig. 7.3.17a shows that even after some growth, the 
delamination is still at the original interface (i.e. between the fourth and fifth 
plies). The view in 7.3.17a is of an area near the y-axis. In Fig.7.3.17b the view is 
of an area offset slightly from the y-axis. There are many ply cracks in this area, 
but the delamination is still between the fourth and fifth plies. 

Fig. 7.3.18 shows a montage of four micrographs of a 90/0/0/90 sublaminate 
which had an initial delamination diameter of 40mm. In this case the delamina- 
tion has grown along the interface between the fourth and fifth plies. After some 
growth, the stresses were such that a crack formed at both ends of the delamina- 
tion. At one end the crack grew all the way to the surface of the laminate. 

The micrographs and the X-rays indicate that even when a configuration starts 
out fairly simple (in this case a single delamination separating the laminate into 
two balanced and symmetric laminates), the damage development tends to be 
complicated. Matrix ply cracking was prevalent in all the tests. The density of 
the cracks varied with location in the sublaminate. The delamination sometimes 
switched interfaces, depending on the stacking sequence. In some cases fiber 
microbuckling occurred. 

The strain at which the delamination began to grow was recorded for each 
configuration, flsing these strains and the analytical results in Figs. 7.3.6-7.3.11 
one can estimate the peak strain-energy release associated with delamination 
growth for each specimen. These results are summarized in Fig. 7.3.19. Both 
G[ and G// are shown. Delamination growth appears to be dominated by G/. 
The critical magnitude of G/ based on Fig. 7.3.19 is not much different than the 
392-513 J/m 2 reported by the materials manufacturer (ref. 46). When G// was 
large, the critical G/ tends to b e so mewhat less. 

For the 30 mm delamination with a (0/90/90/0) stacking sequence, G/ was 
quite low. This was probably due to the high strain which would have been 
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30 mm 40 mm 60 mm 


Fig. 7.3.19 Mode I and Mode II strain-energy release rates for delamination 
growth. Results are shown for two stacking sequences and three delamination 
sizes. 
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required to reach a value of Gj of about 4007/m 2 . Fig. 7.3.20 summarizes the 
strains at which delamination growth began for each configuration. For each 
delamination size the strain was significantly larger for the (0/90/90/0) stacking 
sequence. For the 30mm case the strain for growth was -.0055. This nominal strain 
plus the flexure of the buckled region results in a very high compressive strain in 
some areas. One might guess that this very high compressive strain could trigger 
other damage mechanisms which would augment the effect of the postbuckling 
induced interlaminar stresses. In fact, for both 30mm (0/90/90/0) specimens 
tested, a horizontal crack formed shortly after delamination growth occurred 
(see the X-ray in Fig. 7.3.12). This horizontal crack formed when microbuckling 
occurred on the interior surface of the sublaminate. 

The two stacking sequences differed considerably in terms of the stability of 
delamination growth. The (90/0/0/90) laminate always exhibited at least some 
slow incremental delamination growth with increased load. The delamination 
growth in the (0/90/90/0) laminate tended to be abrupt and often the initial 
growth was quickly followed by formation of a large horizontal crack. 
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Fig. 7.3.20 Strain at which delamination growth began for IM7/8551-7 laminates 
with an embedded delamination. 
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Chapter 8 

SUMMARY AND RECOMMENDATIONS 

There are many mechanisms which contribute to compression failure of lam- 
inated composites. This study has focused on delamination growth driven by 
local buckling of a delaminated group of plies. This mechanism was referred to 
as instability-related delamination growth (IRDG). A review of the literature in- 
dicated that no detailed analysis had been performed for any three-dimensional 
configuration which exhibits IRDG. Detailed stress analysis helps identify which 
material and geometric parameters govern IRDG. Also, without detailed stress 
analyses it is difficult to determine the accuracy of less expensive approximate 
analyses. 

The goal of this study was to enhance the understanding of IRDG in 3D 
configurations through a combination of analyses and experiments. There were 
three tasks in this study: development of a suitable stress analysis program, 
performance of a parametric analytical study, and performance of a combined 
experimental and analytical program. 

A geometrically nonlinear 3D finite element analysis program, NONLIN3D, was 
developed to perform the required stress analysis. This program was tailored to 
exploit the vector processing capability of the CDC VPS32 supercomputer. The 
program uses a 20-node isoparametric hexahedron element. The program was 
designed to perform analysis by substructures. This resulted in a large reduction 
in computational effort, since nonlinear analysis could be restricted easily to just 
that part of a configuration which needed it. An approximate contact analysis 
was also implemented. A variety of checks were performed to assure the accuracy 
of the analysis. 

The parametric analytical study was performed on a fictitious material which 
had the same in-plane stiffnesses as a quasi-isotropic laminate, but had no stacking 
sequence effects. This material was used so that the effects of strain level and 
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geometric parameters could be studied without the additional complications due 
to stacking sequence. Both embedded and edge delaminations were considered. 
These configurations were found to behave in a complex manner. Some of the 
observations and accomplishments were: 

1. The configurations were definitely mixed mode. In some cases G I was larger 
than G[[. For other cases the opposite was true. 

2. In general, there was a very large variation of the strain-energy release 
rates along the delamination front. Usually, one would expect initial delamination 
growth to occur only along a small portion of the front. 

3. The locations of maximum Gj and G u depended on the delamination shape 
and the applied strain. 

4. The mode III component was much smaller than (7/ and G n . 

5. The analysis predicted that for some configurations the sublaminate 
would overlap the base laminate unless contact constraints were added. Most 
of the analyses were performed without imposing contact constraints, since the 
addition of the constraints further complicates an already complicated problem. 
To determine the errors caused by not including contact constraints, a few 
strain-energy release rate results were calculated both with and without contact 
constraints. Tr.ese results showed that for moderate postbuckling strains, the 
errors were fairly small. 

6. The distribution of total strain-energy release rate can be calculated fairly 
accurately using plate analysis. 

7. An intuitive interpretation of the deformation and load transfer for the 
embedded delamination was presented. This interpretation explains much of the 
observed variation of the strain-energy release rates. 

The combined experimental and analytical program involved a variety of 
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configurations. A few preliminary tests were conducted strictly for verifying 
NONLIN3D. The primary goals of this phase of the study was to determine 
whether the behaviors predicted by N0NLIN3D would occur in actual specimens 
and to determine if a strain-energy release rate parameter could be used to predict 
the onset of IRDG. 

IRDG was studied for two material systems: AS4/PEEK and IM7/8551- 7. 
Two sublaminate types were considered: (0/90/90/0) and (90/0/0/90). Because 
of problems in fabricating the PEEK specimens, only qualitative comparisons of 
analysis and experiments were possible. Quantitative comparisons were made for 
the EM7/8551-7 material system. Some of the major conclusions were: 

1. The shape of the postbuckled sublaminate depended on the stacking 
sequence. The predicted shapes agreed well with the experimentally measured 

shapes. 

2. The direction of delamination growth depended on stacking sequence for 
the edge delamination. This behavior was predicted by the analysis. 

3. Delamination growth occurred along only a small portion of the delamina- 
tion front, as predicted by the analysis. 

4. Initial delamination growth in the IM7/8551-7 laminates appeared to be 
governed by the magnitude of G/. 

5. Matrix ply cracking generally accompanied delamination growth. In some 
cases fiber microbuckling also occurred shortly after delamination growth occurred. 

Instability-related delamination growth was defined to be delamination growth 
caused by buckling of a delaminated group of plies. The tests showed that other 
damage mechanisms can also become operative when buckling occurs. Explaining 
the interaction of these various mechanisms with IRDG is beyond the scope of 
this study. In fact, there is so much that needs to be done to fully understand 
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IRDG, that it is difficult to make just a few recommendations for future work. 
The following non-comprehensive comments are offered as general guidelines for 
future work in the area of instability-related delamination growth. 

1. Considerable effort was expended in developing and using NONLIN3D. The 
3D analysis is still quite expensive, in spite of the techniques used to reduce the 
computational effort. These techniques served only to make it possible to make 
the calculations presented herein. If the synergistic effects of various damage 
mechanisms are to be studied with any rigor, new approaches must be developed. 
Hopefully, the results presented will be used to provide direction and insight for 
developing approximate analyses. 

2. A more complete characterization of basic material properties is needed. A 
full 3D set of material moduli is obviously needed for accurate 3D stress analysis. 
Also needed is better interlaminar fracture toughness characterization, particularly 
for mixed mode situations. 

3. The effects of the following on IRDG need to be examined for 3D 
configurations: 

a. unsymmetric sublaminates 

b. global bending 

c. thermal stresses 

d. initial imperfections 

e. multiple delaminations 

4. Synergistic effects due to the presence of multiple damage types need 
considerable experimental investigation. 


185 



REFERENCES 


1. Kachanov, L. M.: Separation Failure of Composite Materials. Polymer 
Mechanics, vol. 12, 1976, pp. 812-815. 

2. Whitcomb, John D. : Analysis of Instability-Related Delamination Growth. 
Proceedings of the Sixth Annual Mechanics of Composites review, AFWAL-TR- 
81-4001, 1981, pp. 95-100. 

3. Whitcomb, John D.: Finite Element Analysis of Instability- Related Delam- 
ination Growth. J. Composite Materials, Vol. 15, Sept. 1981., pp. 403-426 

4. Whitcomb, John D.: Strain-Energy Release Rate Analysis of Cyclic 

Delamination Growth in Compressively Loaded Laminates. NASA TM-84598, 
Jan. 1983. Also, in Effects of Defects in Composite Materials, ASTM STP 836, 
American Society for Testing and Materials, 1984, pp. 175-193. 

5. Whitcomb, John D.: Parametric Analytical Study of Instability-Related 
Delamination Growth. Composites Science and Technology, 25, 1986, pp. 19-48. 

6. Chai, Herzl; Babcock, Charles D. and Knauss, Wolfgang G.: One Dimen- 
sional Modelling of Failure in Laminated Plates by Delamination Buckling. Int. 
J. Solids Structures, Vol. 17, No. 11, 1981, pp. 1069-1083. 

7. Wang, S. S.; Zahlan, N. M.; and Suemasu, H.: Compressive Stability of 
Random Short-Fiber Composites, Part II-Experimental and Analytical Results. 
J. Composite Materials, vol. 19, July 1985, pp. 317-333. 

8. Ashizawa, M.: Fast Interlaminar Fracture of a Compressively Loaded 
Composite Containing a Defect. Presented at the Fifth DOD/NASA Conference 
on Fibrous Composites in Structural Design, New Orleans, LA, Jan. 1981. (also 
available as Douglas Paper No. 6994). 


186 



9. Simitses, G. J.; Sallam, S.; and Yin, W. L.: Effect of Delamination of Axially 
Loaded Homogeneous Laminated Plates. AIAA Journal, 1985, pp. 1437-1444. 

10. Sallam, S. and Simitses, G. J.: Delamination Buckling and Growth of Flat, 
Cross-Ply Laminates. Composite Structures, vol. 4, 1985, pp. 361-381. 

11. Ramkumar, R. L.: Fatigue Degradation in Compressively Loaded Com- 
posite Laminates. NASA CR-165681, April 1981. 

12. Ramkumar, R. L.: Performance of a Quantitative Study of Instability- 
Related Delamination Growth. NASA CR- 166046, March 1983. 

13. Gillespie, John W. and Pipes, R. Byron: Compressive Strength of 

Composite Laminates with Interlaminar Defects. Composite Structures, Vol. 2, 
No. 1, 1984, pp. 49-69. 

14. Konishi, D. Y. and Johnson, W. R.: Fatigue Effects on Delaminations 
and Strength Degradation in Graphite/Epoxy Laminates. Composite Materials. 
Testing and Design, Fifth Conference, ASTM STP 674, American Society for 
Testing and Materials, 1979, pp. 597-619. 

15. Chai, H.: The Growth of Impact Damage in Compressively Loaded 
Laminates. Ph.D. Thesis, California Institute of Technology, March, 1982. 

16. Chai, Herzl and Babcock, Charles D.: Two-Dimensional Modelling of 
Compressive Failure in Delaminated Laminates. J. Composite Materials, Vol. 19, 
Jan. 1985, pp. 67-98. 

17. Webster, John D.: Flaw Criticality of Circular Disbond Defects in 

Compressive Laminates. Center for Composite Materials Report No. 81-03, 1981. 

18. Ashton, J. E.; and Whitney, J. M.: Theory of Laminated Plates. 

Technomic, Inc., Stanford, Conn., 1970, pp. 31-34. 


187 


19. Shivakumar, K. N. and Whitcomb, J. D.: Buckling of a Sublaminate in a 
Quasi-Isotropic Composite Laminate. Journal of Composite Materials, Vol. 19, 
Jan. 1985. 

20. Whitcomb, J. D. and Shivakumar, K. N.: Strain-Energy Release Analysis 
of a Laminate With a Postbuckled Delamination. Fourth International Conference 
on Numerical Methods in Fracture Mechanics. Pineridge Press, 1987, pp. 581-605. 
Also available as NASA TM 89091. 

21. Yin, W. L.: Axisymmetric Buckling and Growth of a Circular Delamina- 
tion in a Compressed Laminate. International Journal of Solids and Structures, 
vol. 21, 1985, pp. 503-514. 

22. Fei, Zhizhong and Yin, Wan-Lee: Postbuckling Growth of a Circular 
Delamination in a Laminate Under Compression and Bending. Proceedings of 
SECTAM XII, Vol. II, Callaway Gardens, May 1984, pp. 130-134. 

23. Bottega, W. J. : The Mechanics of Delamination in Composite Materials. 
Ph.D. Dissertation, Yale University, 1984. 

24. Rhodes, M. D.; Williams, J. G.; and Starnes, J. H. Jr.: Low-Velocity 
Impact Damage in Graphite-Fiber Reinforced Epoxy laminates. Paper presented 
34th Annual Conference Reinforced Plastics/Composites Institute, The Society of 
the Plastics Industry, Inc. New Orleans, La., Jan. 29-Feb. 2, 1979. 

25. Byers, Bruce A.: Behavior of Damaged Graphite/ Epoxy Laminates Under 
Compression Loading. NASA CR-159293, 1980. 

26. Porter, T. R.: Compression and Compression Fatigue Testing of Composite 
Laminates. NASA CR-168023, 1982. 

27. Chai, H.; Knauss, W. G. and Babcock, C. D.: Observation of Damage 
Growth in Compressively Loaded Laminates. Experimental Mechanics, Vol. 23, 
Sept. 1983, pp. 329-337. 


188 



28. Frederick, D. and Chang, T. S.: Continuum Mechanics. Scientific 

Publishers, Inc., Cambridge, 1972, pp. 79-82. 

29. Irons, B. M.: Economical Computer Techniques for Numerically Integrated 
Finite Elements. Int. J. of Numerical Methods in Engineering, vol. 1, 1969, pp. 
201-203. 

30. Pawsey, S. F.; and Clough, R. W. : Improved Numerical Integration of 
Thick Shell Finite Elements. International Journal for Numerical Methods in 
Engineering, vol. 3, 1971, pp. 275-290. 

31. Zienkiewicz, O. C.; Taylor, R. L.; and Too, J. M.: Reduced Integration 
Technique in General Analysis of Plates and Shells. International Journal for 
Numerical Methods in Engineering, vol. 3, 1971, pp. 275-290. 

32. Gentzsch, Wolfgang (editor): Vectorization of Computer Programs with 
Applications to Computational Fluid Dynamics. Friedr. Vieweg k Sohn Verlags- 
gesellschaft mbH, Braunschweig, 1984. 

33. Howser, L. M. and Tennille, G. M.: VPS-32 MINI MANUAL, Central 
Scientific Computing Complex Document S-lb. NASA Langley Research Center, 
1987, pp. 2. 

34. Lambiotte, J.: The Solution of Linear Systems of Equations on a Vector 
Computer. Ph.D. Dissertation, University of Virginia, 1976. 

35. Newman, J. C., JR.: Finite Element Analysis of Crack Propagation - 
Including the Effects of Crack Closure. Ph.D. Thesis, Virginia Polytechnic 
Institute and State University, May, 1974. 

36. Noor, A. K. and Hartley, S. J.: Evaluation of Element Stiffness Matrices 
on CDC STAR-100 Computer. Computers and Structures, vol 9, pp. 151-161. 
1978. 


189 



37. Noor, A. K. and Peters, J. M.: Element Stiffness Computation on CDC 
CYBER 205 Computer. Communications in Applied Numerical Methods, vol. 2, 
1986, pp. 317-328. 

38. Bathe, K.: Finite Element Procedures in Engineering Analysis. Prentice- 
Hall, Inc., Englewood Cliffs, New Jersey, pp. 167-169. 

39. Weaver, W., Jr.: Computer Programs for Structural Analysis, Van 

Nostrand, Princeton, New Jersey, 1967. 

40. Rybicki, E. F. and Kanninen, M. F.: A Finite Element Calculation of Stress 
Intensity Factors by a Modified Crack Closure Integral. Engineering Fracture 
Mechanics, 1977, Vol. 9, pp. 931-938. 

41. Wang, S. S. and Choi, I.: The Mechanics of Delamination in Fibre- 
reinforced Composite Laminates. Part I- Stress Singularities and Solution Struc- 
ture. Part II- Delamination Behavior and Fracture Mechanics Parameters. NASA 
CR 172269 and CR 172270, Nov. 1983. 

42. Raju, I. S. and Newman, J. C., Jr.: Stress-Intensity Factors for a Wide 
Range of Semi-Elliptical Surface Cracks in Finite-Thickness Plates. Engineering 
Fracture Mechanics, Vol. 11, No. 4, p.819. 

43. Timoshenko, S. P. and Woinowski-Krieger, S.: Theory of Plates and Shells. 

0 

McGraw-Hill, New York, 1959, pp. 415-416. 

44. Bostaph, G. M. and Elber, W.: A Fracture Mechanics Analysis for 

Delamination Growth During Impact on Composite Plates. 1983 Advances in 
Aerospace Structures, Materials, and Dynamics, proceedings of the ASME Winter 
Annual Meeting, pp. 133-138. 


190 


C3 



45. Kelkar, A.; Elber, W.; and Raju, I. S.: Large Deflection Behavior 

of Quasi- Isotropic Laminates Under Low Velocity Impact Type Point Loading. 
AIAA/ASME/ASCE/AHS 26th Structures, Structural Dynamics, and Materials 
Conference, 1985. Paper no. AIAA-85-0723. 

46. U. S. /European Aircraft IM7/8551-7 Qualification Test Results. Hercules 
Aerospace Products Group publication, 1987. 


191 



APPENDIX 

SHAPE FUNCTIONS FOR 8- AND 20-NODE ELEMENTS 


This appendix gives the shape functions for the 8- and 20-node elements. The 
node numbering sequence was given earlier in Fig. (2.2.1). 

The shape functions for the 8-node element are : 


5(1) = (1 + <5i) * (1 - S 2 ) * (1 - S 3 )/8 

5(2) = (1 + <$i) * (1 + 62) * (1 — fc)/8 

5(3) = (1 + <5i) * (1 + 62) * (1 + 6 3 )/8 

S( 4 ) = (l+6 l )*(l-6 2 )*(l + 6 3 )/8 
5(5) = (1 - <5i) * (1 — ^2) * (1 — S 3 )/8 

5(6) = (1 - 5i) * (1 + 62) * (1 - «j)/8 

= (1 - ^ l ) * (1 + ^ 2 ) * (1 + 6 3 )/8 
5(8) = (l-*i)*(l- ' ;*(1 + W8 


The shape functions for the 20-node element are : 


Corner nodes: 


5(1) = .125 * (1 - Si) * (1 - 62) * (1 - 6 3 ) * {-Si -62-63-2) 

5(3) = .125 * (1 — 61) * (1 + 62) * (1 — 63) * (— + 62 — S 3 — 2) 

5(5) = .125 * (l — <5j) * (l + 62) * (l + $3) * ( — 4-62 + 63 — 2) 

5(7) = .125 * (l — $i) * (l — 62) * (1 + ^3) * (— — 62 + 63 — 2 ) 

5(13) = .125 * (1 + Si) * (1 - $2) * (1 - S 3 ) * (61 -62-63 - 2 ) 

5(15) = .125 * (1 + Si) * (1 + S 2 ) * (1 - S 3 ) * (6 1 + 62-63-2 ) 

5(17) = .125 * (1 + <$i) * (1 + S 2 ) * (1 + £3) * (Si +62 + 63-2) 

5(19) = .125 * (1 + *0 * (1 - S 2 ) * (1 + 63) * (Si -62+63-2) 

Mid-side nodes on the plane 61 = 0 : 


5(9) = .25 * (1 - Si * Si) * (1 - S 2 ) * (1 - 63) 

5(10) = .25 * (1 - 61 * 61) * (1 + 62) * (1 - *3) 

5(11) = .25 * (1 - * *0 * (1+ 62) * (1 + 63) 

5(12) = .25 * (1 - Si * 61) * (1 - 62) * (1+ S3) 
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Mid-side nodes on the plane 62 = 0 : 


5(2) = .25* (1 - 82 * 6 2 ) * 
5(6) = .25 * (1 — * ^2) * 

5(14) = .25 * (1 — 82 * 62) * 
5(18) = .25 * (1 - 62 * 62) * 

Mid-side nodes on the plane 63 = 0 : 


5(4) = .25 * (1 — £3 * ^3) * 
5(8) = .25 * (1 - S3 * 63) * 
5(16) =.25*(1-63*£ s )> 
5(20) = .25 * (1 - * 3 * *3) 1 


(1 - <50 * (l-« 3 ) 

(1 - ^i) * (1 + 63) 
(1 4- 61) * (1 - $3) 
(1 + 61) * (1 + 63) 


(1 - * (1 + «2) 
(1 - «l) * (l - * 2 ) 

( 1 + 61 ) *(l+fc) 

(1 + «0 * (1 -^2) 
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